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ABSTRACT 


Integral equations (IE’s) are widely utilized to calculate induced currents on anten- 
nas and scatterers, but they are seriously restricted in their ability to handle inhomoge- 
neous penetrable structures having multiwavelength dimensions. The utilization of finite 
element (FE) techniques has not been as pervasive as the use of IE’s. The IE represen- 
tation matrix is “full”, containing few, if any, zero valued elements. The techniques for 
Operating on these large-sized full matrices require undesirable amounts of processor 
time. FE techniques produce sparse matrices due to the strictly local interactions be- 
tween discrete unknowns. The application of FE’s to unbounded problems, however, 
requires supplementary enforcement of the far-field radiation conditions. The Field 
Feedback Formulation (F *) circumvents the full-matrix computational “bottleneck” by 
allowing FE based numerical methods to be employed. Even though the resultant sparse 
matrices may be larger than the “full” matrices discussed earlier, most elements have a 
value of zero. Numerical procedures exist to optimize operations with these sparse ma- 
trices. Calculational speeds can be orders of magnitude faster. Computer techniques to 
implement and validate this new technique are the basis for this thesis. Excellent 


agreement with classical results are demonstrated. 
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I. INTRODUCTION 


A. HISTORY 

The application of finite element techniques to evaluate the solution of differential 
equations is well documented. The utilization of these techniques by _ the 
electromagnetics community has not been as pervasive as the use of integral equations 
(IE's). IE’s are widely utilized to calculate induced currents on antennas and scatterers, 
but they are seriously restricted in their ability to handle inhomogeneous penetrable 
Structures having multiwavelength dimensions. As the complexity and number of nodal 
degrees of freedom grow, the size and dimension of the representative matrix must also 
grow. This matrix is “full”, containing few, if any, zero-valued elements. The available 
numerical techniques for operating on these large-sized full matrices require undesirable 
amounts of processor time. 

Differential equation (DE) based techniques, such as the finite element method. 


produce sparse matrices due to the strictly local interactions between discrete unknowns 





Which result. The application of DE’s to unbounded problems, such as those of scat- 
tering and radiation, require some form of supplementary enforcement of the proper 
far-field conditions. These radiation boundary conditions are innately incorporated into 


integral equations. 


B. FIELD FEEDBACK FORMULATION 

The Field Feedback Formulation (F°) circumvents the full-matrix “bottleneck” in the 
computational process by allowing DE based numerical methods to be emploved. Even 
though the resultant sparse matrices may be larger than the “full” matrices discussed 
earlier, most elements have a value of zero. Numerical procedures exist to optimize 
Operations with these sparse matrices. Calculational speeds can be orders of magnitude 
faster than with full matrices [Ref. 1]. Although the F employs sparse matrices to rep- 
resent the fields in the materials being considered, it does require augmentation to en- 
force the radiation condition at infinity on the scattered fields. This comes in the form 
of a feedback matrix composed of surface integration generated elements. A concept 
evaluation, for a special axisymmetric case, was already accomplished, as detailed in 
[Ref. 2] and [Ref. 3]. Computer techniques to implement and validate this new technique 


are the basis for this thesis. 


C. POTENTIAL BENEFITS 

This thesis will lead to an increased understanding of the advantages and disadvan- 
tages of this novel computational procedure for handling geometrically complex material 
scatterers. This method may ultimately allow computer-aided design of important 
electromagnetic structures such as low-observable aircraft, high efficiency dielectric lens 
antennas, and other electromagnetic scattering occurrences due to atmospheric anoma- 
lies. Structural details and material inhomogeneities, as well as physical dimensions (in 
multiple wavelengths), can be accommodated using the Field Feedback Formulation. 
These capabilities far surpass those which are possible with contemporary integral 


equation techniques for the case of inhomogeneous penetrable scatterers and antennas. 


to 


Il. FORMULATION 


A. INITIAL NOMENCLATURE 
Assume there is a three dimensional object that 1s infinite in one direction. Such an 
object, in cross section could look like Figure 1. This object only varies in two dimen- 


sions, and therefore, is actually a two dimensional (2-D) object. 


V 
TN x 
0S 
Figure 1. A Typical Object 
The wavenumber &, is defined as 
k 2m 27fo 
0 ~~ = eee 

Ao 


Where /, is the free space wavelength associated with an electromagnetic wave of fre- 


quency f, and c is the speed of light. The X and Y Cartesian coordinates are 


wavenumber normalized, such that XY = A,x and Y=4&,y. This coordinate normalization 
will be used throughout this development. A similar technique could also be developed 
in the polar coordinate system. The magnetic field will also be normalized such that 
li = —jy,#, where ny, = = = 120z ~ 377Q, is the impedance of free spacewana 
# is the usual magnetic field in units of A/m. Thus the normalized /7 has the same V/m 


units as E. Potentials may be defined as, 

E,{xy) = ¥(4,¥) (2.1) 
for the transverse magnetic (ITM) case, with £,, £,, and H,=0, and 

Hx y) = W2(4,¥) (2.2) 


for the transverse electric (TE) case, with //,, 7, and E,=0. 

Our objective is to calculate the scattered fields for an arbitrary (2-D) penetrable 
object using the Field Feedback Formulation (F°). As shown in Figure 2, a familiar 
closed loop system illustrates the relationship between the incident, scattered, total and 
far fields. 







Par-Field Green's 





Function Integral 






Oreen's 





Function Integral 
Figure 2. Field Feedback Formulation 


At point 1, the incident field drives the total system. The incident field at 1, combined 
with the scattered field at 4, forms the total field at 2. This total field forms the 
boundary conditions that drive the finite element program at 2. At point 2, initially, the 
incident field drives the U operator. U represents the feed forward operator in the F.. 
This operator uses a finite element technique to solve the boundary value problein. At 


point 3, the boundary conditions are solved for the object perimeter potentials and the 


normal derivative of these potentials. T represents the field feedback operator that takes 
the perimeter potentials and associated derivatives and provides the scattered fields at 
point 4, the offset boundary. The fields at point | and 4 are added (unlike the familiar 
feedback or control system where negative feedback 1s employed). At point 2, a com- 
bined incident and scattered field exists. These fields are the combined boundary con- 
ditions for the finite element boundary value program. These combined fields (total 
fields) are then applied to the L operator to calculate the perimeter fields and the asso- 
ciated derivatives on the objects perimeter. This looping may be repeated until a steady 
State condition, at point 3, 1s reached. The existence of a steady state condition assuines 
stability. Stability for physical svstems should not be a problem. However, when 
mathematically modeled, instabilities may result. The error magnification or condition 
number of the system must also be seriously considered if this iterative looping process 


is to be used. The alternative approach is to form an equivalent svstem, where: 
equivalent operator = U-[I1— T+ U]"', 
with 
f= Identity Matrix 
and 
T-U = Combined Effects of the T and U operators . 
Either approach is viable since. 
ie Winciaenr + T- U + Winciaenr t (T* UY + Winctaens + ©. =LE- Te Ul Winetaent 


This becomes more obvious when w is factored from the equality leaving, 


incident 


AIRS (Uae (GP oO i ae 


Placing this in closed form, 


o, 
5 
i 

a 
| 
rm) | 


Ol — ° rad 
tic to ee 


Finding the equivalent operator will require a matrix inversion and for verv large prob- 
lems this mav lead to excessive computation times. The matrix inversion technique will 
be investigated. The fields may then be extended to anv point in space using a far field 
Green’s function surface contour integral. This far field pattern is available at point 5. 
The incident field 1s usually produced by a plane wave generator. This provides the 
boundary conditions on the offset contour. This contour is called “offset” since it is 
approximately the same “shape” as the objects perimeter but is slightly larger. The dis- 
tance between the perimeter and this contour is called the offset distance and will be 
discussed in Chapter III]. The boundary conditions may be any desired field or wave that 
satisfies Maxwell's equations. These waves may arrive from any direction and be of any 
magnitude. The boundary conditions may also be a composite of any number of waves 
silce superposition does apply to these svstems. A user provided subroutine 1s necessary 
if conditions other than a single plane wave, cvlindrical mode (of arbitrary mode num- 


ber) or individual input boundary condition is desired. 


B. MAXWELL’S EQUATIONS 


Maxwell's equations can be written using our previous normalizations as, 
VxE=n,H 2) 


and 


| 


Vx H=2,E. (2.4) 


“~ 


[Ref, 4 |) Ber o)— — with similar definitions for D,and D,. Equations 2.3 and 2.4 


? 
i 


can be further expanded into the D,, Dy, and Dz components such that. 


u,H, = DyE, (2.5) 
u,H, = —DyE; (2.6) 
u,H, = DyE, — DyE, (2.7) 
¢ E,, = DyH, (2.8) 
¢,E, = —DyH, (2.9) 
¢,E, = DyH, —DyH,. (2.10) 


For the TM case, with propagation in the z direction, 








De. 
H.=— 
and 
Se: 
hor 


Note that the H, field = 0. These two equations can be combined to form, 


v7 ] NN 
H=T-VWy x 2. (2.11) 


Similarly for the TE case, with propagation in the z direction, 


Substituting equations 2.5 and 2.6 into equation 2.10 vields, 





aye 5 a 
Gale. = Dy ie. : (Zaha) 

Substituting equation 2.1] into equation 2.13 yields, 

Dy Dy 4 
E,W, + Dy a, + Dy 7 ==1\); (2.14) 
Equation 2.14 can be further simplified to, 

l Z 
Ve] ZV, | +e =O. (2.15) 


Simularly, equations 2.2, 2.8, 2.9 may be substituted into equation 2.7. This yields, 
I 
Ve] Z-V2 | + ub =0. (2.16) 


Equations 2.15 and 2.16 are TM and TE duals. These two differential equations describe 


the potentials inside the object of interest. Defining _ =o and e,=f for the TM case 


| eae - 
and =-=a and w,=f/ for the TE case and substituting these new definitions into 


equations 2.15 and 2.16 vields one differential equation, 
VeloVy] + fy =0. (2215 


C. VARIATIONAL EQUIVALENCE TO THE DIFFERENTIAL EQUATION 
The Euler-Lagrange variational formulation is based on the stationarity of a func- 


tional, of the function w and its first derivatives. [Ref. 4 ] 


i= | | F(X, Y, u, Vw) dX d¥ (2.18) 


inside § 


It can be shown that the first variation of the functional is zero, 6J=0 , if the 


Lagrangian. F, satisfies the Euler - Lagrange equation: 





C CF C GF Bye 
GAN ( C( Dy) )+ GY ( é(Dyw) aw =I) (2519) 


The problem thus becomes to find the F, which when substituted into equation 2.19, 


vields the original differential equation. The Lagrangian, 
F = o|(Dy) + (Dyb)’} - BY’ 
can be simplified to, 
F=oVw-Vw — Bw’. (2.20) 


When equation 2.20 is substituted into equation 2.19, 





Cr 
= ICBO. 
C(Dyw) xv 
One. 
= De Oe 
C(Dyw) yV 
a 
ai = 2 


Therefore, 


|» 


Qesdya) a (20Dyp) + 2b =0 


B) 


r 
vie | 


oO 


Or, 
D,laDywJ + DylaDyW) + By =0 
which simplifies to, 
V-e(faVy)+ By =0. (2221) 


Therefore, the general functional has been found since equation 2.21 and equation 2.17 


are identical. Substituting the a and f definitions into equation 2.20 yields, 

Fae Vby Vy —e¥7 (2.22) 
and 

Fy =o Vig Vg — lh. (2.23) 


Equations 2.22 and 2.23 are integrated over the interior to S with known boundary 
conditions for either W, or W, on S. To physically interpret the variational formulation. 


it is noted that, 


Vv, =u AHyY - Hy 


i 


and 
Vw = e{EyY a Bat 
aivereiore, 
iG -| | ut Te Baxay 
S 
and 


lee | or »E —p,H «H dX ay. 
S 


Substituting Maxwell's equations into /, gives, 


n=| [ov x E)»H—(V x H)-E dX dy 


S 


n= | [vee x H) dX dY 


Note that E x H is the complex oscillatory Poynting vector. It is different from the 
usual Poynting vector of E x H'. Thus, both functionals, J, and /, are proportional to 
the complex phasors for oscillatory power. The oscillatory power is the phasor repres- 


enting the excess of instantaneous radiated power minus the average radiated power. 


D. FINITE ELEMENT BOUNDARY VALUE SOLUTION 

With the discussion of the variational mathematics completed, a simple rectangular 
boundary value problem will be discussed in detail. The rectangular geometry allows for 
an easier formulation but in no wav limits the solution from being extended to more 


complicated object geometries. Figure 3 on page 11 shows the region of concern. 
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Figure 3. Rectangular Object 


Consider spanning the rectangular region by a triangular mesh, as shown in 
Figure 4 on page 12. Both a and f may be functions of position but may not vary 
within an individual triangular element. Given a fine enough mesh structure, a smooth 
transition in material properties may be approximated. For notational simplicity this 
positional dependance will not be carried forward. The variational approach will yield 
the solution to the boundary value problem by finding the W(x, y) which gives the sta- 
tionary value of, 


a fb 
js (aVy «Vw — By*) dx dy 


which is constrained on the boundary by the previously specified boundary conditions. 
feel, 5] 
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Figure 4. Rectangular Mesh 


The values of w at the interior nodes become the discretized unknowns: 


Wy= WX, Y) f=)... Mandj=1...N 


where 
X,= ie AX 
and 
Yp=jeAY 
for the uniform mesh structure with AY = Trap and AY = —s . Approximating, 
M+1 N41 
W(x, y) = » Wy u(x, y) 
=O.) 


which includes the known boundary nodal values p,, and Wy.,,, for) = O and N + 1, 


and w,, and W, y., fori = Oand M + 1, linear pyramidal basis functions, u,(x, y), Which 


have unit value at the (i,j) node and zero value at all surrounding nodes. Figure 5(a) 1s 
a top view of the pyramidal basis function, while Figure 5(b) is a perspective view. 
1.0 


a. b. 


Figure 5. Pyramidal Basis Function, (a) top view, (b) perspective view 


The functional I will thus be a discrete function of each of the nodal values of y, 
i I(Woo: Wo, toes Wi jp voy Waa, Nei): 


The approximate discrete solution will be found by the system, 


me oan eiiand 1 = 1h. NV. 
Gree 


Now, 
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where, 


M+] N+1 


Vy x, 5) = » » WP wVuj lex, y). 


=O 70 
The gradient of the basis function 1s, 


| avy 
Vein, lS) = By 
m,n 


and the basis function 1s, 


ow 
Urn Aon y= ow 
m,n 





Therefore, the system of equations to solve becomes, 


a fe 
| | (Cone. ~Vy — Bia, nW) ax dy =, for 7t— lee ie 
0 “0 


Substituting for Vwandw in terms of the nodal values of wy _ for 


m=1..,3%andnu=1... i. gives, 


_ 
— 
4 
Pah 


+1 N+] 


d 


aD 
“iy| | CAS n° Vu; = film. ntti, j) dx dy = 0 
ofl 8 lead ¢ 


i=] j=0 
Bene, 

ab 

[ [een © We — Bln nt) do de = Pll, 2), i) 
0 


Regrouping this to put the known nodal values on the right hand side gives for 


m= 1... ¥Pand 7— 14 itenionmmodes: 
MON 
SS yFLon i =- YS WyF Lon, )6)1. 
— Boundary 
Nodes Only 
for Gj) 


14 


By renumbering the nodes using a single index for the unknown nodal values, 
k= NG -—1)+jand/= N0m—-1)+n7 and, 


MeN 


. 


fn 
k=] 


FU, k)by = - > Fl K' Wye. 
rz 


The functional F{/,k) =0 if nodes / and k are not both associated with at least one 
common triangular element. Therefore, Vu,» Vu, and uu, will be zero except in triangles 
Where / and k both appear as nodes. For the example mesh structure of Figure 3 on 
page 11, this produces a banded matrix, F. This sparse matrix can be easily displayed 


by first defining column vectors of the unknown nodal values in the mesh 
Y= Wii Via Mad 
NEL TOMA Ol OCC aN a) ee 


Mies - matrix elements fF [(rn, n).{/,/) pre zero unless the (m.n) node shares at least one 
Bement With the G,}) node. Thus, the node values in w, will be coupled only to 


Ww... W., and any associated boundary nodes. This is written as, 
rae =e [Bae a RGus ans aid [PV 5 


where 4,, Band C,are Vx \ complex arrays and P,isa VN x N, array where .\, is the 


number of boundary nodes associated with the 1-th column. Thus appears as, 


By © Y; P Ys, 
A, B, C, Y> P, ae 
0 Arg) Byey Cuer || Pav ea || See 

Ay By |) Vu ge | ea 


If we denote y, as the initial boundary values and w,,., as the final boundary values 


then ¥, becomes onlv the boundary conditions on the top and bottom at / #0 and 


N+ 1. Note that the above svstem ts tri-block in nature and has a large number of zero 
valued elements. The zero valued elements were omitted for clarity. Each element ts 
actually a matrix and therefore it 1s evident how sparse this system 1s. Each of the 
A,, B.C, and P, matrices are equally sparse, however, a global symmetry does not ap- 


pear. 


E. EVALUATION OF THE F - MATRIX CONTRIBUTIONS 
Given an arbitrarv element as shown in Figure 6, the potential, ¥, can be linearly 


approximated bv [Ref. 5 ], 


Wy 
W(x, y) = (7,57, 1) i (7 ] 4 W2 
V3 
3 
= > betel y), 
L—1 
| 
3 
2 


Figure 6. Mesh Element 


Where w, = W(x,, y,) 18 the nodal value at the k-th node, [7 ]=3 x 3 Transform Array 


and 


ig ~_— 


u(x, y) = Linear basis functions for the k-th node 
Tig 
= (ee i} iG k 


amy . 


13, x 


= les = Ty yy t+ ee 


Note that, 


1 , fork = m 


en 
KX Ym) 0 , fork # im 


It can be shown that, 


(. — ¥3) 3 — V1) (7, — J) 
[7 ]= — (x3 — x) (xy — x3) (x2 — x;) 


(513 — X32) (53h) — 043) G2 — a4) 


where | A] = triangle area, and 


xy Vy I 
2 = cet X4 2 ] 


X3 Ne l 


2A = (xg¥3 + X3V) + XyV2) — (Xgyq + x3 + x2) 


Furthermore, the linear basis functions, u,(x,y), can be interpreted as relative areas of 
the triangle shown in Figure 7 on page 18. 
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Figure 7. Basis Functions Interpretation 


is constant as (x,y) varies and, 


A Ales ¥) 
ts.) = 


Within a given triangular element, the evaluation for k = 1, 2,3 and / = 1, 2, 3 inthe 


element is of interest and, 


| | (,0. © Vu, — Bgtyuty) ax dy. 
inside 
triangle 


q 
The assumption is made that « and £ will be approximated as constants within the tri- 


angle. These material constants can, however, vary from element to element. Taking 


the graditent terms first, 
Vip = Ty fy 


Vu,+ Viz = (Ty x : I hae T> x : T>, ): 


Therefore, 


| | agVu e Vu,dxdy 


triangle 


= 0007) ely t+ Ty, plo, dl Agl, 


where A, 1s the area of the g” element. Next it can be derived that, 


{| ——— A eo Kk ad 
fr | 
utlaxdy = ; ; 
tia , fork=l 


triangle 


More generally, with each yu, raised to an integer power, ,, 


een aT Ny 'No'r3! 
Ui aay ee 
sitar : (m4, + ty + 34+ 2)! 


triangle 


The final result is, 


eG) =| | ory. Vi, — fupt,jdxdyv 


triangle 


q 
Mike 
= AL} aT, aT Uo all) 1? ai jee | 


= AL} a( That T, x) -L 5}. k =. 


F. GREEN’S FUNCTION CONTOUR INTEGRAL 
The scattered fields, w, from an arbitrary object in a vacuum satisfying Helmholtz’s 


equation (see Figure 8 on page 20 ), 
Vt ky =0. 


are {Ref. 6 ], 
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Figure 8. Green’s Function Integration 


where the Green's function 1s, 








and 

a Laman 

a =n+Vw onthe contour 

} 

and 

CG A A JKo (2) 

— e 7 k r= yr’ 
On eee ee ( ol ) 


20 


The Hankel functions of the second kind of order zero and one, He and H) , will 
present a problem for the numeric integration discussed in Chapter V. The imaginary 
portion of these functions rapidly approaches negative infinity as the argument ap- 
proaches zero. The is obtained by a finite difference method using the field 
boundary conditions on surface B (boundary conditions) and the calculated field condi- 


tions on surface P (object perimeter). This results in, 


ow Waconia ~ Wer cier 
offset distance 





p> 
= 


It will be shown that to maximize the accuracy of the Green's function the offset 





| cw 
distance should be made as large as possible. This, however, causes the s to be in- 
Gi 
accurate. Thus, an optimal condition must be found that maximizes the accuracy of the 
entire numeric integration. Such a condition does not maxinuze the accuracy of any one 


of the contributing parts to the Green’s function integrand. 


G. FAR-FIELD EVALUATION 
When the Green’s function integral discussed above is used for far field calculations 
several simphfying relationships develop. These simphfications require a less demanding 


numerical integration. To be in the far field region three conditions must exist, 





l7l>D. I |rl>4, and |¥l>——., 


where /, is the free space wavelength and D is the maximum dimension of the object. 


[Ref. 7 ] As x approaches infinity, 








and 


This requires, 


and 
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= | 
In the far field J+ > ,/-+ andneR—oner. Thus, 





Therefore, 
G(F |r’) as 4) Hot Jor cos @ 
: = 4 thor S 
and 
CGS = Ko 2) —jkgr¢ A A _jkgr’ cos 8 
Ae VARS | ler © (io Fe 


With these new definitions substituted into the original Green’s function integral 


equation, 


= J —jkpr : cw NN ype kyr’ cos@ >, 
Y scattered?) = | Saker OS , Bane + kon er (7 ie : ie. 
RAG 


Cc 


Note that this equation 1s partitioned into a distance dependent term and a theta depend 


term. The theta dependent term mav be defined as 





éy lee 
[= , — +k a) a sr (ee 
Cn 
Cc 
The two dimensional bistatic radar cross section (RCS) per unit length of the cylindrical 
Structure may now be defined as, 


>. ps ? Slee 
RCS = o(¢°, ¢') = lim 22 = | Rises Ae 
TIED ai DUG roo | wy | 


P Wie 








eee er : 
r00 Szhr 4k 


? o ° e e 
where the wavenumber, ky = and the incident field is assumed to be of unit 
‘< “0 

magnitude, | y'| = 1.0. 


Il. MESH GENERATION 


A. INTRODUCTION 

A display or plot of the computer generated mesh structure is not required for the 
problem solution, however, 1t provides an immediate visual confirmation that the in- 
tended problem geometry has been entered correctly into the computer. A large amount 
of initialization data 1s required for even the most rudimentary problem. For this reason, 
the input data is provided to the mesh generation program via a data file. Modifications 


are possible at a later time with minimum effort. 


B. INPUT DATA 

The input data file 1s called INPUT.DAT and contains 25 input fields. Of these 
fields, 14 relate directly to the generation or display of the mesh structure. Only the 
object surface coordinates require adherence to a specific format. All other data need 
only be of the correct type (1.e., character, integer or real). A bref descnptionigigeaem 


field is provided below. 


e Field | is a label of no more than 12 characters. This label is for the plot and input 
data file. 


e Field 21s a character flag that specifies the input coordinate system. If set to “R” 
or “r’, the rectangular system is used. If set to “P” or p . the polar svstemmiseueces 
Ifa“P’."p’.R° or r 1s not detected, an error 1s returned to (hes plae 


* Field 3 is a character flag that if set to “I or 1 will cause several intermeani 
values to be stored to disk during the Finite Element Boundary Value (Reise 
program execution. This option was used during the debug process. 


e Field 4 is a character flag that if set to “D° or “d will create a Din ae 
FORTRAN program capable of replicating the input object. in wavenumber nor- 
malized coordinates. on mainframe computers having a DISSPLA graphics pack- 
age. DISSPLA is a subroutine-based language. The generated program. called 
DISSPLA.FOR, is a compilation of four subroutines calls per element. This file 
can get very large for dense mesh structures. 


* Field 5 1s a character flag that if set to “UL” or “u” will cause a uniform material to 
be assumed. In the uniform case. no material interface exists. This option was 
only used to verify the Finite Element solution accuracy. 


e Field 6 is a character flag that if set to “M” or “m” will cause only the mesh to be 
generated. This is very useful when first starting a problem and the optimal mesh 
structure has not been determined. 


e Field 7 is a real number specifving the desired mesh resolution in wavelengths. The 
mesh resolution determines the dimension of the mesh elements. 


Field 8 is a real number specifving the distance, in wavelengths, between the object 
perimeter and the offset boundary contour. 


Field 9 is a real number that specifies a multiplicative scaling factor for the nu- 
merical integration stepping function discussed in Chapter V. 


Field 10 is a bias term that can be used to shift the numerical integration stepping 
function. This term is used for distances less than 1.0. 


Field 11 1s a bias term that can be used to shift the numerical integration stepping 
function. This term is used for distances greater than 1.0. 


Field 12 is a real number specifving the maximum distance bevond which no further 
contribution to the Green’s Function Integral is made. If this feature is not de- 
sired, this term) should be made larger than the objects maximum dimension plus 
twice the offset distance. 


Field 13 is an integer specifving the number of input data points. 


Field 14 is an integer specifving the angular resolution, in degrees, desired for the 
final radar cross section calculation. 


Field 15 1s an integer specifving the mesh generation technique. 


Field 16 is an integer specifving the perimeter node from which the bisection seg- 
ment originates. This node is called the “start node”. 


Field 17 is an integer specifving the perimeter node on which the bisection segment 
terminates. This node is called the “stop node”. 


Field 18 is a pair of real numbers (on two lines) specifying the x and v coordinates 
by which the object will be displaced. 


Field 19 is a real number specifving. in wavelengths, the desired distance between 
the origin and the first input data point. Fields 18 and 19 when used together, al- 
low an object to be placed at anv position and scaled to anv size. 


Field 20 is a pair of real numbers (on two lines) specifving the real and imaginarv 


parts of for the ITM case or 7 fOmthe hE case. 


Field 21 is a pair of real numbers (on two lines) specifving the real and imaginarv 
parts of ¢, for the TM case or u, for the TE case. 


Field 22 is a character flag that if set to “P” or “p” enables a plane wave generator. 
The plane wave 1s propagating down the vy axis, and generates an £,e’” condition 
on the offset boundary contour. If the flag is set to “C” or “c”, a cylindrical mode 
generator 1s enabled. This generates an E, cos m@ condition on the offset boundary 
contour. If a “P” , “p”, “C” or “c” is not detected, then manually input boundarv 
conditions must follow, and fields 23 and 24 are not used. 


Field 23 1s a real number specifying the wave amplitude. 
Field 24a is a real number specifying the wave frequency (in Hz). 
Field 24b (only for cylindrical case) is an integer specifving the mode number, n. 


Field 25 is the object perimeter data, in either polar or rectangular form. 


An example of an input data file for a homogeneous circular cylinder is provided in 
Appendix A. The majority of the input data is echoed to the computer display and a 
system “pause” is initiated to allow for user inspection. The program may be aborted 
or continued at this ttme. The initial object dimensionalization has already occurred, 
and the number of unknowns and the maximum unknown width is displayed. These 
factors give an excellent indication of the expected run time for the FEBV routines. For 
example, a problem with 512 unknowns and a maximum unknown width of 31 took 806 
seconds to execute while a run with $ unknowns and a maximum unknown width of 3 
took only 10 seconds to execute. These times are for a Intel 80386 based personal 
computer with an Intel 80287 co-processor chip calculating the fields for a circular cyl- 


inder. The FEBV routines are the next code block to execute after the pause is cleared. 


C. MESH GENERATION PROGRAM 

The mesh generation program consists of seven subroutines. These subroutines are 
an integral part of the finite element program and, therefore, were not separated. These 
routines are discussed below. 

l. YO (Input/Output) 

This subroutine reads the information contained in the INPUT.DAT file dis- 
cussed earlier. A two dimensional object can be described in any number of ways, 
however, for simplicity, the polar and rectangular coordinate systems are used. In either 
case, the initial assumption is that all data points are referenced to a local origin. This 
local origin can be offset bv any desired amount using field 18. This offset is independ- 
ent of the entered data points or any size scaling provided by field 19. A plot label file, 
named TEXT.LBL, is created and the initial object perimeter (coded for a display in 
blue) is written to the output file, PLT.DAT. These data points describe the perimeter 
of the object. This will later prove helpful in determining the conformity of the gener- 
ated mesh to the input perimeter. All subsequent screen writes are coded for a display 
in green. Two example objects are shown in Figure 9 on page 27. These objects will 
be used throughout this chapter. The circular cylinder was generated by a separate 
computer program. The “horseshoe” shaped object was manually input. Graph paper 
was used to determine the x and v coordinates of the 28 unequally spaced perimeter 


points. 
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Figure 9. Typical Objects 


2. Rotate 

This subroutine reorders the input data points to allow for any desired bisection 
seginent start and,or stop node. This subroutine ts only used if manual selection of the 
bisection segment start and/or stop nodes 1s requested. 

3. Bound (Boundary) 

This subroutine sub-divides the object perimeter based on the mesh resolution 
specified 1n Field 7. The mesh resolution is the approximate length, specified in wave- 
lengths, that the user desires the perimeter to be divided into. The division of the per- 
imeter is based on linear interpolation between input data points. A new perimeter node 
is placed at each of the sub-division points. Additional input data points are necessary 
in areas of rapid change to allow for a correct object perimeter representation. Ihe ob- 
ject bisection extends from the “start node” to the “stop node”, These nodes are user 
specified in Fields 16 and 17 and, in general, divide the object in half. The bisection 
should be arranged such that the object width, perpendicular to the bisection segment, 
is miminuzed. Thus, a long slender object should be oriented for a major axis bisection. 
Whether the major axis is ortented vertically or horizontally is of no importance. This 


is demonstrated in Figure 10 on page 28 (right side). 
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Figure 10. Bisection Segment and Row Arrangement 


For more complex objects, the bisection is not a straight line, but rather a series 
of lmne segments that approximate a curve. The number of perimeter nodes on the Icft 
side of the bisection segment must equal the number of perimeter nodes on the right side 
of the bisection segment. The bisection segment also contains this same number of 
nodes. Thus, the program inust adjust the requested mesh resolution to ensure proper 
nodal spacing. Thus allows for a piecewise continuous segment to cross the object. 

4. Normal 

This subroutine calculates the unit normals to the object perimeter at each 
newly established perimeter node. The normal ts a perpendicular constructed to the 
chord connecting the two nodes adjacent to the node for which the calculation 1s oc- 
curring. This perpendicular originates at the current node and is of unit length. See 


Figure 11 on page 29. 
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Figure 1!. Unit Normal Calculation 


5. Nodset (Node Set) 

This subroutine uses a two sweep technique to compute the number of nodes 
on each nodal row. The rows are labeled 1 = 1, 2, 3, .... | maximum (IMX) consec- 
utively from the top of the object to the bottom. Segments normal to the perimeter 
surface, calculated in NORMAL, connect the perimeter nodes to the offset boundary 
contour. Iwo such completed normal segments, one on each end, complete a nodal row. 
When all rows are complete, a second contour has been established that approximates 
the object’s perimeter. This contour is displaced from the perimeter by the distance 
specified in Field 8, This contour provides the boundary conditions for the FEBV rou- 
tines. The optimal selection of this offset distance will be a topic of Chapter IV. 

With the object divided into rows, each row can be further divided into equally 
spaced nodes. Based on the requested resolution and whether the objects dimensions 
are expanding or contracting, the nodal spacing is adjusted to keep the elements 
approximatelv the same size. Two adjacent nodal rows form an elemental row. These 
rows are given the same label, ] = 1, 2, 3, ..., IMX as the upper nodal row. The nodes 
On an elemental row are numbered consecutively, starting with the left-most node. This 


node is always an offset boundary contour node. When the upper nodal row is 
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numbered, the process 1s continued for the lower nodal row. An example of this element 
row numbering scheme is shown in Figure 12 on page 30. 


8 16 
15 


13 14 


Figure 12. Element Row Numbering Scheme 


A mesh orientation attribute is set for the left and right portion of each element row. 


This attribute determines if the object has a clockwise or counter-clockwise orientation 


See Figure 13 on page 31. Switching between the four available mesh orientations al- 


lows for a mesh that more accurately conforms to the input object perimeter without 
resulting in a disproportionate mesh structure. 
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Figure 13. Mesh Orientation Attributes 


There is an additional requirement that the number of nodes on a given left or right half 
row must be only one more, equal to, or one less than the number of nodes on the ad- 
jacent left or night half row. Thus, a two sweep process is utilized. This process ensures 
that this requirement is met and that the first and last row have only two nodes. The 
second and second from last row (if present) must have three nodes. It is possible, as 
shown in Figure 14 on page 32, to generate a mesh that has only three rows. This object 
Was Input as a circular cylinder. It is evident that, due to a small number of elements, 
the generated mesh more accurately resembles a square. This will be a problem for cal- 
culating the scattered fields, however, the internal fields can be accurately approximated. 
In this special case, the second and second from last rows are the same. It will be shown, 
in Chapter V, that since this mesh does not closely approximate the input objects per- 
imeter the resulting Green’s function contour integrals and subsequent far field calcu- 
lations have reduced accuracy. 

The nodes of the nodal rows forin the vertices of triangular elements. An ele- 
ment, as defined in Chapter II, has three vertices, each assigned a unique number (1, 2 
or 3). The ordering of these vertices depends on the mesh orientation attribute. Each 


element also has an associated relative dielectric constant (¢€,) and relative permeability 


(4t,). 
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Figure 14. Three Row Cylinder 


6. Sorter 
This subroutine generates a complete mesh row and all the element/node inter- 
connection relationships. With the nodal structure in place, individual nodes can be 
assigned to individual elements within the global mesh structure. Each element in a 
given row ts assigned a unique local element number starting with the left most element 


(relative to the bisection segment). See Figure 15. 
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Figure 15. Element Numbering Scheme 
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It is vital that a method of determining which nodes form the vertices of a given element 
and which elements have a vertex attached to a given local node. The nodes that are 
not part of the offset contour have unknown field values. A single node can be con- 
nected to as many as six elements, only four of which can be in the current row. This 
relationship is shown in Figure 16(a). This hexagonal arrange:nent of six elements 1s 
very common within the mesh. Along the bisection segment or where the mesh orien- 
tation attribute changes from one row to another, this pattern 1s disrupted. An example 
of an extreme case (the center of a circular cylinder) is shown in Figure 16(b). After all 
elements and nodes in an elemental row are assigned, an ordered sweep is conducted of 
that elemental row to determine which nodes are connected to which elements, which 
elements are connected to which nodes, and how manv elements a single node is con- 


nected too. A complete row has now been generated. See Figure 17 on page 34. 
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Figure 16. Two Possible Element Intersections, (a) extreme, (b) normal 


The information necessary to generate the rows and elements will be used again, in 
Chapter IV, to solve the (FEBY) problem. 
7. Finder 
This subroutine determines the x, y coordinates of each node. This data is also 
necessary to solve the FEBV program of Chapter IV and provides the plotting 
coordinates for the PLT.DAT file. This process is repeated for all rows. Thus, the entire 
object 1s generated as the compilation of elemental rows made of triangles. See 


Figure 18 on page 34. 
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Figure 17. Typical Row Structure 
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Figure 18. Global Mesh Structure 


A file 1s now available for display using a commercially available program called 
“CURVE-DIGITIZER”. Any program that can accept x, y coordinate data will 
accomplish the same display process. Minor changes to the MESH program provided 
in Appendix B may be necessary since several “"CURVE-DIGITIZER” plotting codes 


are embedded in the file generation code. 
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D. OPTIMIZATION OF THE MESH 

For most objects, 1t is recommended that the MESH program be executed in the 
mesh generation onlv mode (Field 6 set to “M” or “m”) prior to the execution of the total 
finite element program. This will ensure that the desired mesh structure 1s obtained prior 
to solving for the unknown field values. The MESH (generation only) program requires 
onlv a few seconds for even the most dense mesh Structures. It is readily seen that as 
the mesh resolution increases, the number of rows increases linearly while the number 
of unknowns increases geometrically. Therefore, the mesh calculation times may not 
change appreciably when the mesh is made more dense, however, the calculation time 
for the finite element program will rise geometrically. For this and other reasons, the 
mesh density should be kept low. This will lead to a smaller number of unknowns and 
result in faster program execution times. Figure 19 on page 36 plots this relationship 
for a circular cvlinder. The solution for these unknown field values will be a topic of 
Chapter IV. 

Six different mesh generation methods are available to create mesh structures. Some 
of these methods were evolutionary in nature and provide limited practical benefit. 
Methods | and 6 are by far the most useful. Each method 1s discussed below. 

Method | constructs the bisection segment by connecting the first input data point 
to the midpoint of the perimeter. This segment is divided into equal length segments 
each separated by a node. This method 1s useful for very simple objects such as the 
circular cylinder shown in Figure 20 on page 37. This circular cvlinder will be used in 
the explanations of all mesh generation methods. These illustrations are not intended 
to optimize the mesh structure but rather allow for easy comparison of the 6 basic 
methods. 

Method 2 constructs the bisection segment as described in method 1. The bisection 
segment nodes are, however, ordered differently. This segment is divided by connecting 
line segments from corresponding nodes on the left and right side perimeter. Where 
these line segments intersect the bisection segment, a node is placed. This leads to un- 


equally spaced bisection segments. This can be seen in Figure 21] on page 37. 
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Figure 19. Unknowns and Number of Rows Versus Mesh Resolution 


Method 3 modifies method | by allowing the user to specify, using Field 17, the stop 
node for the bisection segment. This method can be useful to rapidly adjust around 
slightly irregular objects. This can be seen in Figure 22 on page 38. 

Method 4 combines methods 2 and 3 by using the connected line segment technique 
of method 2 to determine the node positions on the bisection segment, the stop node of 
Which is specified as in method 3. This can be seen in Figure 23 on page 39. 

Method 5 improves upon method 4 by repositioning the unequally spaced bisection 
nodes. Linearly interpolated positions for the nodes leads to equal spacing. This 
method reduces the node “bunching” that frequently occurs with methods 3 and 4. This 


can be seen in Figure 24 on page 40. 
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Figure 20. Method ! Mesh Structure Example 
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Figure 21. Method 2 Mesh Structure Example 
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Figure 22. Method 3 Mesh Structure Example 


Method 6 is the final improvement in which method 5 was modified to allow for a 
user-specified start node. This provides the user with the abilitv to start and stop the 
bisection segment at any input node without having to rearrange the input data. Since 
method 6 contains all of the capabilities of the other five methods, tt 1s almost exclu- 
sively used for mesh generation. This can be seen in Figure 25 on page 40. There are 
situations that could be best served by one of the other methods. For example, a cir- 
cular cylinder and other simple symmetric objects can be represented using method lI. 
Method 2 would work well with a square or diamond shape. These objects would be 
bisected by a diagonal. Method 3 would be most suited for a symmetric object with a 
planar material interface. The more dense mesh would be used for the higher 
pernuttivity material. Method 4 would work well with a rhombus or other shghtly 
asymmetric object. Method 5 is almost as versatile as method 6, and would work well 
in any situation where the start node is fixed. A great deal of planning is not needed in 
designing most mesh Structures since they are calculated and displaved in a matter of 


seconds. 
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Figure 23. Method 4 Mesh Structure Example 


Iterative selection of different mesh generation parameters has proven to be the best 
technique. A summary of all mesh generation capabilities is provided in Table 1. 
Appendix C contains a program called READ.FOR. This program takes the output 
data file from the “Curve-Digitizer” CAD program, called FINALDWG.DAT, and after 
receiving the answers to several prompted questions, creates a new INPUT.DAT file. 
Typically, the CAD program is used to generate the perimeter of the object. This may 
be as a series of points, line segments or a combination of the two. The answers to the 
prompted questions provide the additional information needed to fill the remaining data 
fields. This is intended to be a first step towards allowing a user to specify or design an 
object and then be able to calculate the scattered fields from this object. Although it is 
far from efficient, it does serve a definite purpose. Further refinement of this program 
will allow for an easier user interface. This should enable technically trained personnel, 


without a detailed understanding of these programs, to benefit from the Field Feedback 
Formulation. 
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Method 5 Mesh Structure Example 


Figure 24. 
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Method 6 Mesh Structure Example 


Figure 25. 
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Table 1. SUMMARY OF MESH GENERATION CAPABILITIES 
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IV. FINITE ELEMENT BOUNDARY VALUE PROGRAM 


A. INTRODUCTION 

The Finite Element Boundarv Value (FEBV) program is the feed forward (U) oper- 
ator in the Field Feedback Formulation (F) as shown in Figure 2 on page 4. This 
program solves the Helmholtz equation, as discussed in Chapter II, for the Dirichlet 
boundarv condition specified in the input data file, as discussed in Chapter III. These 
boundary conditions are imposed on the offset boundarv contour. The purpose of this 
program is to find the unknown field values inside the offset boundary contour within 
the input object. Since the ultimate goal is to obtain the scattered far fields from the 
object, the unknown field values on the perimeter are of primary interest. It will also 
be necessary to approximate the normal derivative of the field at the object perimeter. 


As derived in Chapter II, the goal of this program is to solve, 
LAI, + LBAY; + LGigae == Balper 


The A matrix represents the effect that the I] — 1 row has on the 7” row values. Simularly, 
the B matrix represents the effect that the I-th row has on itself, while the C matrix re- 
presents the effect that the I + | row has on the I-th row. All other rows do nomened 
the I-th row. This 1s due to the pyramidal basis functions, discussed in Chapter I], that 
all have zero value when a node is not directly connected to another node. The P matrix 
represents the combined effects of the boundary nodes on the I-th row. These values 
are transposed to the right side of the equality to form the system forcing function. It 
is Worth remembering that for the I = 1 row, the A matrix = 0 and that for the I = 
IMX row, the C matrix = 0. Thus, using the row by row stepping process, first dis- 
cussed in Chapter III of the mesh generation process, the A, B, C and P matrices can 
be filled. There is an A, B, C and P matrix for each row, and it is necessary to calculate 
the functionals for two elemental rows to fill one set of matrices. It is, therefore, obvious 
that the data is used twice, once for the current row and again for the next row. Spe- 
cifically, the functionals necessary to complete the B matrix and totallv fill the C matrix 
is used again to totally fill the next row’s A matrix and partially fill the B matrix. Thus, 
fora row I - 1, (such that 1 - 1 4 1), all of the A matrix and part of the B and P matrices 
are filled. When the row is stepped to I, the B and P matrices are completed and all of 


the C matrix is filled. The A, B, C and P matrices represent tri-block matrices within a 


much larger global matrix. This row of tri-block matrices in the global matrix can be 
partially solved using the forward portion of the Ricatti transform. The Ricatti trans- 
form is a numerical technique that optimizes the solution for banded (tri-block) systems 


of linear equations. The equations are, 
Ri) = — (B+ AR) C, 
and 


Sie) = (8B, = (B+ AiR) ' « (P — A,S)) 


t 


Sameres 1s the 7*+ 1 R matrix and the S_, is the “+ 1S vector. Note that the next 
rows R matrix and S vector depends on the previous rows R matrix and S vector. 
During the forward step (MARCH subroutine) an R matrix and S vector is generated 
for each row. When all rows have been calculated, a back sweep computes the unknown 


field values, y,.. The equation is 
Wi = Rib; + Sp 


This back sweep must read the RS data from the disk backwards. This is a storage in- 
tensive process for all but the most trivial problem. The field values wW,_, is first found 
Meeemicmberine that the last row (1 = IMX ). C,,,=0. With Cryy=9 , Kins.) = 0 
Gy — Sry,7-, . Which 1s the last calculated S$ vector. The recursion continues until 


all of the wy, ‘s have been calculated. 


B. FINITE ELEMENT BOUNDARY VALUE PROGRAM 
The eight subroutines comprising the FEBV program are discussed below. 
1. Zero 
This subroutine fills all the arrav positions of the A, B, C and P matrices with 
QO + j0. This zeroing is necessary after all calculations concerning a given row are 
completed. 
2. Varint (Variational Integration) 
This subroutine calculates the complex functionals for a given input element. 
The functionals reflect the effect that each node has on the three nodes associated with 
an element. The subroutine must be provided with the X, Y coordinates of the element 
nodes, and the material parameters ¢, and y,. These functionals are returned in a 3 x 3 


complex matrix. The area of the input element is also provided as a by-product of this 
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numerical integration. The areas of each element are sorted and used to find the largest 
and smallest elements. Once found, the largest and smallest areas are combined to form 


the area ratio, which 1s defined as, 


maximum area 


area ratlo = —— 
minimum area 


This ratio should be kept as small as possible. For a circular cylinder, values slightly less 
than 2.0 are possible. For more complicated objects, the ratio can get considerably 
larger. A ratio => 2.5 causes a screen message of YOU SHOULD CONSIPa 
ABORTING THIS RUN AND LOOKING AT THE MESH IN CURVE DIGITIZER. 
A BETTER METHOD MAY BE AVAILABLE”. It may or may not be possible to 
obtain an area ratio < 2.5. The intention of the area ratio is to insure that a uniform 
mesh is constructed prior to attempting a problem solution. Smaller area ratios are in- 
dicative of mesh structures that do not have grossly different element sizes. This will 
lead to more accurate finite element solution. 
3. Fill 

This subroutine calculates and stores all of the functionals for an entire ele- 
mental row. This 1s accomplished by repeated VARINT calls for each element that 
comprises an elemental row. 

4. BNDC (Boundary Condition) 

This subroutine calculates and stores the boundary conditions desired for the 
offset boundary contour. These conditions can be plane wave, cylindrical modes or in- 
put manually. For plane wave conditions, the percent error of the FEBV program will 
also be returned. This calculation onlv has meaning for the uniform case (Field 5 set to 
“CL” or “u’). Memory limitations do not allow this feature for the cylindrical mode 
boundary conditions for modes other than zero and one. A Separate routine could be 
made available for offline percent error calculations. No such feature is provided if the 
boundary conditions are manually input. Any desired boundary condition may be gen- 
erated by modifving or appending the proper code to this subroutine. 

5. Loader 

This subroutine loads the A, B, C and P matrices for a given row. For all but 
the 1 = 1 and IMX rows, two calls of FILL/VARINT for two elemental rows of data 
are required to fill the A, B. C and P matrices. This is an additive process that starts 
after the ZERO subroutine initializes the matrices. Successive functionals are added to 


the existing data in the proper array positions. 


a. 


— a —$<antty = 


— 


6. March 

This subroutine performs the forward portion of the Riccati transform. The 

output data, called RS.DAT (R matrix and S vector), are stored on disk. 
7. CSMINV (Complex Square Matrix Inversion) 

This subroutine accomplishes the complex matrix inversion required by the 
forward Riccati transform. A maximum dimension of 50 x 50 was established to limit 
memory utilization. This allows for a maximum unknown width of S50. 

8. Sweep 

All of the above subroutines are called at least once for each row. The sweep 
routine is called only after all of the row calculations are completed, and then onlv once. 
This subroutine conducts the Riccati back sweep by reading the RS data generated by 
the MARCH subroutine. This data is read off the disk backwards, using the 
FORTRAN “backspace” command. The returned field values are actually individual 
contributions due to a unit valued basis function being individually applied to each of 
the boundary nodes. In so doing, the problem need only be solved once. After the data 
is stored. in matrix form as the L.DAT file, the boundarv value problem may be solved 


for any incident field by multiplying the U matrix by the new incident fields. Thus. 
Vinernees =[U]- Vi incident’ 


A summarv of the input, output and error data is provided in the OUTPUT.DAT file. 
9. Save 

This subroutine was necessary to allow for reasonably sized problems. Since 
all of the Field Feedback Formulation code could not fit into 640 kilobytes of memory, 
the program was divided into two parts. All of the data necessarv to perform the pro- 
grams is saved in the F3.DAT file. This data is then read bv the field feedback program. 
This technique, though verv inefficient, circumvents the 640 kilobyte memory limitation 
imposed by the IBM Disk Operating System (DOS). Efforts to convert this code to run 
under a compiler that does not have a 640 kilobyte memory limitation, such as Micro- 
way NDP FORTRAN compiler, will be pursued at a later date. 


C. VARINT VALIDATION 

The first step in the program validation required an understanding of the error 
convergence of the variational elements as a function of element size, d and material 
properties, ¢,andu,. The test program provided in Appendix D varies the element size, 


in wavelengths, for a test mesh structure. This test structure shown in Figure 26 on 


page 46 has only one unknown at the center. The mesh is made of four adjacent ele- 
ments. These elements are inside a uniform material having properties, ¢, and y,. The 


other four nodes are established as boundary nodes. 


PLANEWAVE 





Figure 26. Test Mesh Structure 


A plane wave, of user specified frequency and amplitude, is used to determine the 
boundary conditions on the four boundary nodes. This plane wave is propagating down 
the vertically oriented axis. The unknown nodal value is calculated and compared to the 
actual plane wave value (I + j0). A percent error is calculated as the dimension of the 
elements, d are reduced. As expected, the solution became more accurate as the ele- 
ments become smaller. A plot of the results of this test program is provided in 


Figure 27 on page 47. This data was generated with ¢,= 1+ 0, andy, = 1+,0. 
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Figure 27. Solution Error for a Test Mesh Structure 


The region of interest is where the error approaches zero. An expanded plot of this re- 
gion is provided in Figure 28 on page 48. The accuracy of the solutions were desired 
within | percent of the exact value. This requires an element that is < ~ , in the ma- 
terial. To ensure the desired error was achieved, elements were usually scaled to be 
< a , in the material. The phase of the plane wave was also varied to determine the 


im 20 
effect on convergence. No significant effect was noted. 
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Figure 28. Solution Error for a Test Mesh Structure (Expanded) 


D. FINITE ELEMENT BOUNDARY VALUE PROGRAM VALIDATION 

The next validation step required an actual object to calculate the fields in and on. 
The simplest object was actually not an object at all, but rather an imaginary circular 
cvlinder in free space. A plane wave was propagated through this free space. Since no 
material interface existed, the exact solution would be known for all positions. The 
plane wave established the boundary conditions on the offset contour. Trial runs were 
conducted varying the mesh resolution and boundary offset contour distance. With the 


error defined as, 





2 
X(calculated value — actual value) 
ero. = | ——— 


E(actual value)’ 


two errors were defined. The perimeter error only considers the perimeter nodes 1n the 
error calculation. The bisection segment error only considers the bisection segment 


nodes, with the exception of the two end nodes, in the error calculation. The end nodes 
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of the bisection segment are part of the perimeter and are, therefore, not considered. 
As expected, the perimeter error could be rapidly reduced by decreasing the offset dis- 
tance. This is due to the fact that the node or nodes closest to an unknown node dom- 
inate the field contributions at this node. Again, a goal of < 1% error was desired. The 
reduced offset distance had a very small effect on the bisection segment error. The only 
Way to significantly reduce this error was to increase the mesh resolution. An increased 
mesh resolution reduced both errors. Figure 29 shows the perimeter error as a function 


of the number of bisection segment nodes. 
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Figure 29. Perimeter Error 


The three curves are for offset distances of 0.05, 0.025 and 0.01 4. Figure 30 on page 
50 shows the bisection segment error as a function of the number of bisection segment 
nodes. The two curves are for offset distances of 0.05 and 0.025 /,. For offset distances 
less than 0.025 A, , the curves are almost identical to the 0.025 2, curve. These curves 


are onutted for clarity. 
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Figure 30.  Bisection Segment Error 


The calculated and exact field values for a 0.5 A, diameter circular cylinder, with a 
0.054, mesh resolution and offset distance and ¢, = 1 +0 is shown in Figure 31 on page 
St. The exact fields values for the real and imaginary portion of the plane wave are 
shown as squares and diamonds, respectively. The solid curves are the calculated field 
values. The perimeter and bisection errors were 0.74 and 1.69 percent, respectively. 
Figure 32 on page 52 shows the effects of not properly adjusting the mesh resolution 
and offset distance when the material is changed. In this case, the permittivity was 
changed to ¢,=4+/ 0. This caused the perimeter and bisection segment errors to 1In- 
crease to 18.2 and 41.1 percent, respectively. The mesh resolution and offset distance 
were reduced to 0.0217, and the permittivity was changed to ¢,=4—j4. The results are 
shown in Figure 33 on page 53. This proper selection of mesh resolution and offset 
distance reduced the perimeter error and bisection segment error to 0.44 and 0.56 per- 


cent, respectively, Note that in the lossy material case the two errors are very close in 
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magnitude. This is due to the way that the error is defined. This definition, in a lossy 


material, overemphasizes the objects leading edge. 
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Figure 31. 0.5 / cylinder, c, = 1 + j0 


This portion of the object (as well as the trailing edge) is the most accurate for the 
bisection segment calculations. This is because of the relative closeness of these nodes 
to the perimeter. For the lossless material, the bisection segment error is tvpically two 
to three times the perimeter error, for the same number of bisection nodes. The con- 
clusion that the offset distance should be made as small as possible in order to minimize 
the perimeter error is correct, but will prove to be counterproductive in the long run. 
There is a competing effect that will require this contour offset distance to be as large 


as possible. This effect, associated with the accuracy of the Green’s function contour 


integral, will be discussed in Chapter V. 
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Figure 32. 0.5 / cylinder, ¢, = 4 + j0 
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Figure 33. 0.5 1 cylinder, ¢, = 4 — (4 


E. INHOMOGENEITY 

Until this point, all testing was done in circular homogeneous dielectric materials. 
It was, at a minimum, desirable to place the object in a vacuum to calculate the scattered 
fields. This requires the first two and last two elements of each row to have 
€é-=1+4+/,/0 and vp,=1+/0. These elements represent the space between the objects per- 
imeter and the offset boundary contour. Since the plane wave solution is no longer valid 
for this geometry, a new problem with a known solution was needed. 

Given a homogeneous dielectric circular cylinder of radius a and permittivity ¢, , as 


shown in Figure 34 on page 54, it can be shown that [Ref. 8 ], 


WAR, 6) = A,J,(4,R) cos nd, 
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Figure 34. Cylindrical Mode Geometry 


Monere, 


w,(R, @) is the field value inside the dielectric material. 


Sal Bye 


n= 
nav ,, 


TAR In(kpRa) 1 (Rg) — kyl” (ky Rg) HR,)] — 


_ HOR Ink, Ra)F' (Ra) ~ KJ’ (kK Ra)I (Ra) J 
J, = Bessel Function of the 1° kind of order n 


and 
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An additional formulation is required for the fields outside of the cylinder. A cylindrical 
wave (mode, n = 1!) was applied to a homogeneous dielectric cylinder. The exact and 
finite element field values were calculated and compared. Figure 35 on page 55, shows 
a typical result. Mesh resolution, offset distance and material permittivity were varicd 
to determine convergence. The results were similar to the dielectric homogeneous plane 


Wave case discussed earlier. 
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Figure 35. Typical Cylindrical Node Boundary Value Problem Result 


F. FINITE ELEMENT CONCLUSIONS 

High accuracy solutions can be obtained by maintaining a mesh resolution of 
<< in the material. The offset distance should be kept at approximately the same 
magnitude as the mesh resolution, but may be as large as a Increased accuracy re- 
quires an increased mesh resolution. This increase in accuracy is at the expense of in- 


creased processor time. 
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V. FIELD FEEDBACK PROGRANI 


A. INTRODUCTION 
The Field Feedback Program 1s the feedback (T) operator in the F . This program 
utilizes the output of the finite element program and the input incident fields to evaluate 


a “near field” Green's function integral. This integral is called “near field” in that the 





integral will be performed within of the object perimeter. As seen in Figure 36, 


i) 
20 = a ° e 
three surface contours are defined. The boundary contour 1s where the incident field 
boundary condition is applied. 


eo Geometric Contour, GC 
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Figure 36. Contour Arrangement 


The perimeter contour is where the finite element program solves for the field values on 
the object perimeter. The geometric contour is midway between the boundary and the 
perimeter and is the contour over which the Green’s function contour integral is per- 


eae Cl : ; _ 
formed. This 1s necessary to allow the ah to be approximated by a finite difference 
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Hecnamique. [he w,- 1s the average of the field values on the boundary and perimeter. 
Thus, 


Ow (wy boundary ar 





On On 
and 
(y perimerer + V boundary) 
nt i> — 


The Field Feedback program 1s provided as Appendix E. 


B. FIELD FEEDBACK PROGRAM 
1. Input 
This subroutine reads the finite element data stored in the F3.DAT file by the 
SAVE subroutine. All necessary nodal X. Y coordinates are calculated. 
2. IMAT (T Matrix) 
This subroutine loads a complex matrix called TMAT. Each element of the T 
matrix is the result of a Green’s function contour integral. This integral is conducted 
on the GC contour. The m-th column of the T matrix is evaluated with a single unit 


valued basis function on the m-th boundary node. The equation, 


esc 


may now be numerically evaluated at each of the n boundary nodes. This process 1s 








CW a 
votes On a 
GCC 


repeated for each row until the entire matrix 1s filled. 

The numeric integration uses a rectangular mid-point approximation technique. 
For this linearized problem, this is equivalent to a trapezoidal integral technique. The 
integral path between any given two nodes is subdivided based on the distance to the 
first point. This subdivision mavbe controlled by three input variable fields. A multi- 
plicative scale factor and offset terms are available. To date, the utilization of the scale 
factor (set > 1.0) and the offset terms (set > 0 ) have not proved necessary. This may 


be necessary for more complicated geometric structures. 
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3. CNSOLV 


This subroutine solves the equation, 
1 
Ce =([]—- heigl r Wy ek 


where C, is the total scattered fields on the boundary, I is the identity matrix and 
Wooundery 1S the initially specified incident fields. 
4. FFLD (Far Fields) 
This subroutine calculates the scattered far field or bistatic radar cross section 
per unit length of the object. Any desired integer angular resolution for the calculation 
may be specified. The final results are stored in the FFPAT.DAT file. 


C. FIELD FEEDBACK VALIDATION 
Given a plane wave boundary conditions imposed on an imaginary cylinder in free 
Space, a Green’s function contour integral may be performed around the cylinder. This 
cvlinder is imaginary in the sense that it does not actually exist. The cylinder is actually 
an artificial boundarv that does not form a material interface. For points inside this 
cvlinder, the resulting integral should equal the negative of the actual plane wave value 
at that point in freespace. For points outside the cylinder, the resulting integral should 
equal zero. Several test conditions were investigated to ensure that the numerical inte- 
gration of the Green’s function did arrive at the expected values. Maximum absolute 
errors for these tests were less than 7 x 10". These test cases started With exact values 
of the field. w and the normal derivative, oS . After initial validation, finite element 
calculated w and oa values were used. Errors increased by approximately a factor of 
10. Numerous competing effects were observed with two dominant effects becoming 
evident. 
I. Small Object Phenomenon 
For very small objects, where only a few elements are needed to meet the = 
mesh resolution requirement, the normal derivative accuracy 1s crucial. Since the normal 
derivative is the calculated normal to the piecewise linear approximation of the objects 
perimeter, this fit is usually not adequate. To prevent this problem, additional 
perimeter'boundary nodes are needed. This is easily accomplished by increasing the 
mesh density, which is controlled by the mesh resolution (input field 7). For these very 
small objects, the RCS is almost uniform for the TM case and co-sinusoidal for the TE 


case. The utility of this calculation must, therefore, be questioned. 
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2. Offset Distance Phenomenon 
The offset distance must not be made too small or the accuracy of the Green's 
function integral will suffer. For a circular cvlinder, the optimal offset was between 0.03 


; jae / 
4 and 0.035 4. Note that these distances are always < ay 


VI. VALIDATION 


A. INTRODUCTION 

The difficulty in the total Field Feedback Formulation (F) program validation 1s in 
finding problems that have well established or accepted solutions. The total F * program 
is the combination of the mesh generation’finite element program and the field feedback 


program. If possible, a problem that has a closed form analytical solution 1s desired. 


B. HOMOGENEOUS CIRCULAR CYLINDRICAL SCATTERING 

This is the first test case for the F° program. A penetrable circular cylinder has a 
closed form analytic solution, so the final results could be verified against exact values. 
It can be shown that for the TM case (E - wave) the reflection coefficient 1s, 
yb Rad y(Re)— af GET alba Re) 
pre _ 


n = 


5 Ge 
lb RMA (Ry) — FE Jl by Red) 


and that for the TE case (H - wave) the refection coefficient is, 


Ln(b RoW (R,) = [FE Ilka (R a 
p™ 


TR” (R= Je — Sk RHR) 


— [ €, ; 
where A. = /e., },=./> and —= /— . (Ref. 9 ] The scattered field fommire 
IM case 1s, 


EAR, 6) = — gies ry + 2) iT HOR) cos ro 
n=] 


Similarly, the scattered field for the TE case 1s, 


H3(R, 6) =— eon! rat + 2) iT HOR) cos ro 


n=] 
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In the far field region, as R approaches infinity, the scattering width may be defined as, 


| | 77°|? 
o(¢) - 2n lim. e | pyincident | 2 


and 


4 OO 
Ce) ies IT) +2) T, cos nd |?, 


n=] 


The source code, without the Bessel function routines, for these calculations is provided 
as Appendix F. 

Three different radius dielectric cylinders were tested. Wavenumber normalized radi 
oO and 2.0 with ¢, = 2.56 + j0 were selected. The TE and TM results of the 0.5, 
1.0 and 2.0 radu problems are provided as Figure 37 on page 62, Figure 38 on page 
63, and Figure 39 on page 64. To validate the Chapter VI conclusion that the unit 
normal was the cause of the errors for very small circular cvlinders, the ¢, was increased 
from 2.56 to 25.6. This value still meets the requirement for mesh resolution not to ex- 
ceed a0 . The TE and IM results are provided as Figure 40 on page 65. Actual cross 
section values are indicated bv the squares and the Field Feedback Formulation sol- 


utions are plotted as a single solid curve. 


C. HOMOGENEOUS IRREGULAR OBJECTS 

The results detailed in [Ref. 10] were next validated with the F program. This ob- 
feemeds seen in Figure 41 on page 66, is a dielectric shell with inner radius = .25 A). 
outer radius = .30 /, and e, = 4 + j0. The comparison of the [Ref. 10] data and the 


3 ‘ ; = 
fF results are provided as Figure 42 on page 67. 
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Figure 37. Cylinder, TE and TM Case, k,a = 0.5, €, = 2.56 
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Figure 38. Cylinder, TE and TM Case, k,a = 1.0, &, = 2.56 
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Figure 39. Cylinder, TE and TM Case, k,a = 2.0, &, = 2.56 
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Figure 40. Cylinder, TE and TM Case, k,a = 0.5, ¢, = 25.6 
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Figure 41. Dielectric Shell Mesh 


The TE case shows excellent agreement. The ITM case tends to diverge at the smaller 


values of @. 


D. AN INHOMOGENEOUS OBJECT 

The results detailed in (Ref. 11] were next validated with the F program. This ob- 
ject, as seen in Figure 43 on page 68, is a dielectric ring with inner radius = .25 /, , outer 
radius = .30 A, and e,= 4 + j0. The exact solution to this problem is available. [Ref 
11] 
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Figure 43. Dielectric Ring 


The mesh generation program was not modified to conform to both the outer and inner 
material boundaries. The original mesh generation program design concept was to 
closely fit the mesh to the objects outer perimeter and then allow for material inhomo- 
geneities to be accounted for by different material parameters being assigned to individ- 
ual elements. A section of the generated mesh is shown in Figure 44 on page 69. Note 
the additional curved line showing the .25 7, inner radius. This curve does not follow 
the established element boundaries. The effective ring is shown in Figure 45 on page 
69. This resulted in the inner radius having an irregular pattern as elements near the 
material interface varied in material composition. This is similar to the granular noise 
problem characteristic of delta modulation communication channels. The TE and 1M 


results are provided as Figure 46 on page 70. 
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Figure 44. 


Figure 45. 
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Figure 46. Dielectric Ring, TE and TM Case, ¢, = 4 + j0 
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VII. CONCLUSIONS 


A. RESULTS 


The Field Feedback Formulation proved to be an excellent tool for calculating the 


internal and scattered fields from the tested inhomogeneous asymmetric objects. The 


kevs to satisfactory results are, 


, 
Keep the maximum element dimension < 30° 
Maintain a mesh with a uniform distribution, and 
Maintain the offset distance near the mesh resolution in magnitude, but no greater 


9 
than x0 


The techniques used proved to be capable of adapting to a wide varietv of situations. 


These adaptations did not require program modifications or reprogramming. The nu- 


merous built-in features, as discussed in the input field description of Chapter III, al- 


lowed for this robustness. 


B. RECOMMENDATIONS AND EXTENSIONS 


With the basic program validated, future testing and development should, 


Emphasize large object validation against accepted results. 
Stress inhomogeneity and irregularity in all testing. 


Optimize the T matrix element calculation by improving the Green’s function 
contour integral. 


Evaluate the usefulness of the maximum distance feature (Field 12). Bevond this 
miaxinium distance no contribution is made to the Green’s function integral. It is 
not expected that this will be possible for objects less than a few wavelengths in 
maximum dimension. 


Modify the mesh generation program to allow for conformity to multiple interfaces 
(multi-lavered objects without the granular noise problem). 


Modifv all programs for Intel 80386 based Fortran compiler use, such as NDP 
FORTRAN by Microway. This will remove the 640 KB memory limitation im- 
posed by the IBM DOS and allows for the solution to larger scattering problems. 
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APPENDIX A. INPUT DATA FILE EXAMPLE 


CIRCLE = PLOT iaReEr 
R - RECTANGULAR INPUT DATA 
NI - NO INTERMEDIATE DATA RECORDING 
ND - NO DISSPLA PROGRAM GENERATION 
NU - NON UNIFORM MATERIAL (MATERIAL INTERFACE PRESENT) 
NM - NO MESH GENERATION ONLY 
D019 - REQUESTED MESH RESOLUTION (CIN WAVELENGTHS ) 
O03 - CONTOUR DISTANCE (IN WAVELENGTHS ) 
118, - MULTIPLICATIVE GFI SCALE FACTOR 
0 - DISTANCE < 1.0 BIAS TERM 
0 - DISTANCE > 1.0 BIAS TERM 
OD 92 - MAX DIST BEYOND WHICH THERE IS NO CONTRIBUTION TO GFI 
36 - NUMBER OF INPUT DATA POINTS 
i2 - NUMBER OF POINTS FOR SIGMA CALCULATION (CIN THE CIRCLE) 
- ELEMENT GENERATION METHOD 
- START NODE 
5 = s1OP NODE 


= XA AALS TORPSES 

= Y AAI S OEE SES 

- DESIRED DIMENSION (ORIGIN TO FIRST POINT ,WAVELENGTHS ) 
REAL PART OF ALPHA 

- IMAGINARY PART OF ALPHA 

- REAL PART OF BETA 

- IMAGINARY PART OF BETA 

- PLANEWAVE GENERATOR ENABLED 


Wr WIORORrR OF FW Ue 
OG) i.e). 
1 


. 00000000 
. 08682409 
LPO LOG 
. 25000000 
toclouso0 
PoosU2220 
. 43301270 
. 46984630 
. 49240390 


. 25000000 = 


- PLANEWAVE AMPLITUDE 
- INPUT FREQUENCY OF THE PLANEWAVE, FO 
. 50000000 
- 49240390 
- 46984630 
. 43301270 
; DO5UzZ220 
CVE Selo ral, 
eZo000000 
. 17010 
. 08682407 
. 50000000 = 
. 49240390 = 
. 46984630 = 
. 43301270 =, 
~o0 202220 = 
532159560 ar 
. 25000000 = 
Ae EF OE OG, 3) ss 
. 08682404 =i, 
. 00000004 ae 
. 08682413 =| 
ly OT 10 =" 
- 43301270 
132159860 =. 
, 902024220 =. 


00000002 
08682411 
Py cba 
25000000 
32139360 
38302220 
43301270 
46984630 
49240390 
50000000 
49240390 
46984630 


Souec7 20 
GL G)s) Gren) 


- ANGLE, X, Y COORD. 


743501270 
. 46984630 
. 49240390 
. 90000000 
. 49240390 
. 46984630 
. 43301270 
Poo 2022 20) 
3239500 
. 24999990 
le LOLOOO 
. 08682401 


= 25000000 
=.1/ 101000 
-. 08682403 
. 00000007 
. 08682416 
lO 1010 
ne U000T0 
2159390 
5 SUA 
. 43301280 
. 46984630 
- 49240390 
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APPENDIX B. READ PROGRAM 


REAL A(0: 2000), B(0:2000), MRES, DIST, XORIGIN, YORIGIN 

REAL C, D, E, F, DPER, DD, FREQ, EO, MAXD 

INTEGER I, J, METHOD, STARTND, STOPND, NRES, MODE, LBIAS, GBIAS 
CHARACTER*1 CHAR,CHAR1,CHAR2 ,CHAR3,CHAR4 ,CHARS 

CHARACTER"13 NAME 


OPEN(CUNIT = 1, FILE 
OPEN(C UNIT = 2, FILE 
J=0 


'D: FINALDWG. DAT' ) 
'C: MSFORT INPUT. DAT' ) 


WRITE(*,**) ‘ENTER THE PLOT NAME OR LABEL (MAX OF 13 CHARACTERS) ' 
READ(*,1005) NAME 

WRITE(*,*) ‘ENTER THE COORDINATE SYSTEM IN USE. (R OR P)' 
READ(**, 1000) CHAR 

WRITE(*,*) ‘DO YOU WANT INTERMEDIATE VALUES ? (I) ' 
READ(**,1000) CHAR1 

WRITE(*,**) ‘DO YOU WANT A DISSPLA PROGRAM GENERATED ? (D)' 
READ(*,1000) CHAR2 

WRITE(*,*) ‘UNIFORM SLAB (NO MATERIAL TO VACUUM INTERFACE ? (U)' 
READ(**, 1000) CHAR3 

WRITE(*,**) 'MESH GENERATION ONLY ? (M)' 

READ(**, 1000) CHAR4 

WRITE(**,**) ‘ENTER THE MESH RESOLUTION. (0.4, ETC...) /' 
READ(*,**) MRES 

WRITE(*,**) ‘ENTER THE CONTOUR DISTANCE. (0.3, ETC... )' 
READ(**,**) DIST 

WRITE(*,*) 'ENTER THE GFI SCALE FACTOR. (1.0, ETC...)' 
READ(*,**) DPER 

WRITE(*,*) 'ENTER THE GFI ( < 1 BIAS TERM (0,1, ETC...)' 
READ(*,**) LBIAS 

WRITE(*,**) 'ENTER THE GFI ( > 1 BIAS TERM (0,1, ETC...)' 
READ(**,**) GBIAS 

WRITE(*,**) 'ENTER THE GFI MAX DISTANCE (999.0, ETC...)' 
READ(**,*) MAXD 

WRITE(*,*) 'ENTER THE NUMBER OF POINTS FOR SIGMA CALC. (36,ETC... )' 
READ(*,**) NRES 

WRITE(*,%*) ‘ENTER THE DRAWING METHOD. (1 - 6)' 

READ(**,**) METHOD 

WRITE(*,*) 'ENTER THE STARTING NODE NUMBER. (MUST BE > 0) ' 
READ(**,*) STARTND 

WRITE(*,*) 'ENTER THE BISECTION STOPPING NODE NUMBER. ' 
READ(**,**) STOPND 

WRITE(*,*) 'DESIRED DISTANCE FROM ORIGIN TO POINT 1 (IN WLs) ' 
READ(**,*) DD 

WRITE(*,**) ‘ENTER THE X AXIS ORIGIN. ' 

READ(*,**) XORIGIN 

WRITE(*,*) ‘ENTER THE Y AXIS ORIGIN. ' 

READ(*,*) YORIGIN 

WRITE(**,*) 'ENTER THE REAL COMPONENT OF ALPHA. ' 

READ(* ,**) e 
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WRITE(*,*) 'ENTER THE IMAGINARY COMPONENT OF ALPHA. ' 
READ(*,*) D 
WRITE(*,*) ‘ENTER THE REAL COMPONENT OF BETA. ' 
READ(*,*) E 
WRITE(*,*) ‘ENTER THE IMAGINARY COMPONENT OF BETA. ' 
READ(* ,**) F 
WRITE(*,*) ‘PLANE WAVE (P) OR CYLINDRICAL MODE (C) ' 
READ(**,1000) CHARS 
WRITE(*,*) "ENTER THE WAVE FREQUENCY (IN HERTZ) ' 
READ(* ,*) FREQ 
WRITE(*,*) ‘ENTER THE WAVE AMPLITUDE ' 
READ(*,*) EO 
IF( (CHARS. EQ. 'C'). OR. (CHARS. EQ. 'c')) THEN 
WRITE(*,*) 'ENTER THE MODE NUMBER ' 
READ(*,**) MODE 
ENDIF 
mOr10, I = 0, 1999 
READ(1,*,ERR = 20) ACJ), BC(J) 
IF(A( J). GT. 1000) THEN 


colo ss 
Pehl oletOO0)= THEN 
GOTO 5 
ELSE 
J= J 1 
ENDIF 
CONTINUE 
CONTINUE 
a- J - 1 


WRITE(2,1005) NAME 
WRITE(2,1000) CHAR 
WRITE(2,1000) CHAR1 
WRITE(2,1000) CHAR2 
WRITE(2,1000) CHAR3 
WRITE(2,1000) CHAR4 
WRITE(2,1010) MRES 
WRITE(2,1010) DIST 
WRITE(2,1010) DPER 
WRITE(2,1020) LBIAS 
WRITE(2,1020) GBIAS 
WRITE(2,1040) MAXD 
WRITE(2,1020) J 
WRITE(2,1020) NRES 
WRITE(2,1020) METHOD 
WRITE(2,1020) STARTND 
WRITE(2,1020) STOPND 
WRITE(2,1010) XORIGIN 
WRITE(2,1010) YORIGIN 
WRITE(2,1010) DD 
WRITE(2,1010) C 
WRITE(2,1010) D 
WRITE(2,1010) E 
WRITE(2,1010) F 
WRITE(2,1000) CHARS 
WRITE(2,1010) EO 
iaGCGHARS. EQ. © ). OR. (CHARS. EQ. 'c’)) THEN 
WRITE(2,1020) MODE 


LS 


ENDIF 
WRITE(2,%*) FREQ 
DO" 30. = ee) 
WRITE(2,1030)0 1 ACL) BU 
CONTINUE 
WRITE(2,*) 'END OF FILE’ 


FORMAT(A1) 

FORMAT(A13) 

FORMAT(F8. 6) 

FORMAT(I3) 

FORMAT( 1X, 13,7X,F12. 8,5X,F12. 8) 
FORMAT(F8. 3) 


CLOSE( 1) 
CLOSE( 2) 
STOP 

END 
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APPENDIX C. MESH GENERATION/FINITE ELEMENT PROGRAMNI 


FINITE ELEMENT MESH GENERATION PROGRAM 
WRITTEN BY LT. T.B. WELCH 


w/ PROGRAMMING IDEAS FROM PROF. M.A. MORGAN 
COMPLETELY FILE C INPUT. DAT) DRIVEN 
INPUT PARAMETERS: 


CHARACTER (POLAR OR RECTANGULAR INPUT DATA - FLAG) 
CHARACTER C INTERMEDIATE MATRICES WRITE - FLAG) 
CHARACTER (DISSPLA PROGRAM GENERATION - FLAG) 
CHARACTER CUNIFORM MATERIAL - FLAG) 
CHARACTER (MESH GENERATION ONLY - FLAG) 
MESH RESOLUTION 
CONTOUR OFFSET 
DESIRED SCALE FACTOR FOR GREEN'S FUNCTION INTEGRAL 
Pao eConle nh) rOR GFI STEP WHEN < 1.0 
BIAS (SHIFT) FOR GFI STEP WHEN > 1.0 
MAXIMUM DISTANCE BEYOND WHICH THERE WILL BE NO 
CONTRIBUTION TO THE GREEN'S FUNCTION INTEGRAL 
NUMBER OF POINTS 
NUrteneon POINTS FORSCRUSS SECTION, CIN THE CIRCLE) 
MESH GENERATION TECHNIQUE 
START NODE 
OTOP NODE 
X ORIGIN 
Y ORIGIN 
DESIRED DISTANCE (FROM ORIGIN TO FIRST POINT, WAVELENGTH) 
ALPHA (INPUT AS A AND B THEN CONVERTED TO COMPLEX ALPHA) 
BETA (INPUT AS A AND B THEN CONVERTED TO COMPLEX BETA) 
CHARACTER (PLANEWAVE GENERATOR - FLAG) 
AMPLITUDE 
FREQUENCY CIN HERTZ) 
=wORy=< 
AMPLITUDE 
CYLINDRICAL MODE NUMBER (GENERATOR ALWAYS DOES N=1) 
FREQUENCY CIN HERTZ) 
so OR =< 
BOUNDARY CONDITIONS 
DATA POINT PAIRS (X, Y OR RADIUS, THETA) 
END OF FILE MARKER 


COMPLEX ALPHA,BETA,BCOND( 100) , LINE(50) ,ANS( 100) 

COMPLEX A(50,50),B(50,50) ,C(50,50) ,P(50, 100) , SURBC( 100) 

COMPLEX U( 100,100) ,PSI(50,100) ,SVEC(50, 100) 

REAL MRES,MESH(0: 200,5),NDP(200,2) ,OFFSET,DPER ,MRESW 

REAL MINAREA, MAXAREA ,AREA,EO,XORIGIN, YORIGIN, KO ,MAXD 

INTEGER NPNTS,PERND ,NODES( 200) ,MOR( 200) , BIND,NBMX,UNK, LBIAS,GBIAS 


TE 


INTEGER LND(0: 200,3) ,NDL(200,4) ,NCT( 200) ,MAXEL, IMX ,MODE ,NRES 
INTEGER METHOD ,STARTND,STOPND , MINROW , MINEL, MAXROW , NABC( 100, 3) 
CHARACTER*1 CHAR, CHAR1, CHAR2, CHAR3, CHAR4, CHARS 

EQUI VALENCEMGES1 , F) 


COMMON/BLK1/MESH 

COMMON/BLK2/PERND, BIND 

COMMON/BLK3/NODES 

COMMON/BLK4/LND ,NDL,NCT 

COMMON/BLK5/NDP 

COMMON/BLK6/MOR 

COMMON/BLK7 /MINAREA , MINROW , MINEL, MAXAREA , MAXROW , MAXEL, AREA 
COMMON/BLK8/A,B,C,P 

COMMON/BLK9/CHAR2, CHAR3, CHAR4 


OPEN (UNIT = 1, FILE = ‘MATRIX. DAT’ ,STATUS = 'UNKNOWN' ) 
OPEN (UNIT = 2, FILE = ‘PLT. DAT' ,STATUS = 'UNKNOWN') 
OPEN (UNIT = 3, FILE = INPUD DAD soh{lUss— ) Orpem) 

OPEN (UNIT = 4, FILE = 'TEXT. LBL' ,STATUS = 'UNKNOWN') 
OPEN (UNIT = 7, FILE = 'RS.DAT’,STATUS = 'UNKNOWN’ , 


CACCESS=' SEQUENTIAL’ , FORM='UNFORMATTED' ) 

OPEN (UNIT = 8, FILE = PSI. DAT' , STATUS = 'UNKNOWN' , 
CACCESS= ' SEQUENTIAL’ » FORM=" FORMATTED’ ) 

OPEN (UNIT = 9, FILE = 'ABCP. DAT’ , STATUS = 'UNKNOWN' ) 


OPEN (UNIT = 10, FILE = 'BCOND. DAT', STATUS = " UNKNOWN’ ) 
OPEN (UNIT = 11, FILE = 'PMATRIX. DAT’ , STATUS=' UNKNOWN’ ) 
OPEN (UNIT = 12, FILE = 'NABC. DAT’ ,STATUS=' UNKNOWN’ ) 
OPEN (UNIT = 13, FILE = ‘OUTPUT. DAT’ ,STATUS=' UNKNOWN’ ) 
OPEN (UNIT = 19, FILE = 'ABRMAT. DAT’ ,STATUS='UNKNOWN' ) 
OPEN (UNIT = 20, FILE = 'DISPLA. DAT’ ,STATUS=' UNKNOWN’ ) 
OPEN (UNIT = 30, FILE = 'U. DAT‘ ,STATUS=' UNKNOWN’ ) 

OPEN (UNIT = 40, FILE = 'F3.DAT' ,STATUS=' UNKNOWN’ ) 
MINAREA = 999999, 9 


MAXAREA = -999999.9 


CALL IO(NPNTS ,MRES ,METHOD, STARTND,STOPND, OFFSET, ALPHA, BETA, BCOND, 
CEO ,CHAR1,MODE, XORIGIN, YORIGIN, CHARS , DPER, KO, NRES, MRESW, LBIAS,GBIAS 
C ,MAXD) 

CALL ROTATE(NPNTS , METHOD , STARTND , STOPND) 

CALL BOUND(MRES ,NPNTS , METHOD , STOPND ) 

CALL NORMAL 

CALL NODSET( METHOD ,NABC , NBMX , UNK) 

PAUSE 'PLEASE PRESS ENTER TO CONTINUE, OR CTRL BREAK TO ABORT! ' 

CALL LOADER( BCOND , OFFSET, ALPHA, BETA ,NABC, IMX,NBMX,E0,SURBC,CHAR1, 
CLINE ,MODE , XORIGIN, YORIGIN, SVEC, CHARS) 

CALL SWEEP( IMX,NABC,SURBC, LINE ,CHAR1,U,BCOND, PSI, ANS ,CHARS) 

CALL SAVE( BCOND,ANS,U,OFFSET, PERND, CHARS ,DPER,KO,XORIGIN, YORIGIN, 
CNRES , MRESW, LBIAS, GBIAS ,MAXD) 


CLOSE(1) 
CLOSE(2) 
CLOSE(3) 
CLOSE(4) 
CLOSE(7) 
CLOSE(8) 
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CLOSE(9) 

CLOSE( 10) 
MOSE( 11) 
CLOSE( 12) 
CLOSE( 13) 
CLOSE( 19) 
CLOSE( 20) 
CLOSE( 30) 
CLOSE(40) 


STOP 
END 


SUBROUTINE IOCNPNTS ,MRES ,METHOD,STARTND,STOPND , DIST ,ALPHA, BETA, 
CBCOND ,EO,CHAR1 , MODE, XORIGIN, YORIGIN, CHARS , DPER, KO, NRES , MRESW, 
CLBIAS ,GBIAS ,MAXD) 


THIS SUBROUTINE READS IN THE INPUT PARAMETERS 
AND SURFACE DATA POINTS. THESE POINTS CAN BE 
IN EITHER POLAR OR RECTANGULAR FORM. 


COMPLEX ALPHA,BETA,BCOND( 100) 

REAL A,B,PI,XORIGIN, YORIGIN,E0,K0,FO,SFAC,DD 

REAL MESH(0: 200,5),MRES,THETA, RADIUS , DIST,C, LAMBDA, MRESW 

REAL DISTW,DPER,MAXD 

INTEGER I,II,J,K,NPNTS,RES,METHOD,STARTND ,STOPND , MODE, NRES, LBIAS 
INTEGER GBIAS 

CHARACTER*1 CHAR,CHAR1,CHAR2,CHAR3,CHAR4, CHARS 

CHARACTER*12 NAME 


COMMON/BLK1/MESH 
COMMON/BLK9/CHAR2, CHAR3, CHAR4 


C = 2.997925E+08 
PI = 4. O*ATAN(1. 0) 
READ(3,1070) NAME 
READ(3,1000) CHAR 
READ(3,1000) CHAR2 
READ(3,1000) CHAR3 
READ(3,1000) CHAR4 
READ(3,1000) CHARS 
READ(3,%*) MRESW 
READ(3,%*) DISTW 
READ(3,%*) DPER 
READ(3,%) LBIAS 
READ(3,%*) GBIAS 
READ(3,°*) MAXD 
READ(3,**) RES 
READ(3,**) NRES 
READ(3,**) METHOD 
READ(3,*) STARTND 
READ(3,%) STOPND 
READ(3,**) XORIGIN 
READ(3,*) YORIGIN 
READ(3,*) DD 
READ(3,*) A 
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READ(3 57") 9B 

ALPHA = CMPLX(A,B) 
READ(3,*) A 
READ(3,*) B 

BETA = CMPLX(A,B) 
READ(3,1000) CHARI 


IF(METHOD. EQ. 1) THEN 
WRITE(6,*) ‘METHOD 1 
ELSEIF( METHOD. EQ. 2) THEN 
WRITE(6,*) 'METHOD 2 
ELSEIF(METHOD. EQ. 3) THEN 
WRITE(6,*) 'METHOD 3 
ELSEIF( METHOD. EQ. 4) THEN 
WRITE(6,*) 'METHOD 4 
ELSEIF( METHOD. EQ. 5) THEN 
WRITE(6,*) 'METHOD 5 


SELEC IEG 
SELECTEDS= 
SELECTED - 
SELECTED - 


SELECTED - 
ELSE 
"METHOD 6 


WRITE CO.) SELECTEDS- 


ENDIF 


IF( (CHAR2. EQ. 1°) OR. (CHARZ BO): ia 
WRITE(6,*) ‘INTERMEDIATE MATRIX 
ELSE 
WRITE(6,*) ‘INTERMEDIATE MATRIX 
ENDIF 


IF((CHAR3. EQ. 'D'). OR. (CHAR3. EQ. 'd')) 
WRITE(6,*) 
ELSE 
WRITE(6,*) 
ENDIF 


IF((CHAR4. EQ. 'U'). OR. (CHARG. EQ. 'u')) 


OBJECT BISECTION 


CONNECTED NODES 


SELECTED STOP NODE 


t 


CONNECTED/SELECTED STOP NODE’ 


EQUALLY SPACED CONNECTD NODE’ 


SELECTED START AND STOP NODE' 


THEN 
FILE GENERATION 


FILE GENERATION 


THEN 


"DISSPLA FORTRAN PROGRAM GENERATION 


"DISSPLA FORTRAN PROGRAM GENERATION 


THEN 


ENABLED' 


DISABLED ' 


ENABLED ' 


DISABLED' 


WRITE(6,*) 'UNIFORM MATERIAL SPECIFIED (NO INTERFACE) ' 


ELSE 
WRITE(6,*) 
ENDIF 


IF( (CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm')) 
WRITE(6,*) 'MESH GENERATION <<< 
ELSE 
WRITE(6,*) ‘MESH GENERATION AND 


ENDIF 


IF((CHARI-EO.) P ). OR. (GHARIT- EOuap op 
READ(3,*) EO 
READ(3,*) FO 
LAMBDA = C/FO 
KO = 2*PI/LAMBDA 
WRITE(6,**) 


THEN 
ONY 22 


FE PROGRAM 


THEN 


"PLANEWAVE BOUNDARY VALUE GENERATION 


"MATERIAL SPECIFIED WITH A VACUUM AROUND OBJECT' 


- ENABLED' 


- ENABLED ' 


- ENABLED ' 


WRITE(6,*) ‘AMPLITUDE(EO) = ', EO,', WAVENUMBER(KO) = ‘ ,KO 


WRITE(6,*) ' WAVELENGTH — 


" ,LAMBDA,', FREQUENCY(FO) 


ELSEIF( (CHARI. EQ. 'C'). OR. (CHAR1. EQ. 'c')) THEN 


READ(3,**) EO 
READ(3,**) MODE 


80 


= ho 


ep Yai i a in 


10 


20 


READ(3,*) FO 
LAMBDA = C/FO 
KO = 2*PI**LAMBDA 
WRITE(6,*) 'CLYINRICAL BOUNDARY VALUE GENERATION - ENABLED’ 
WRITE(6,*) 'AMPLITUDE(EO) = ', EO,’, MODE NUMBER = ', MODE 
WRITE(6,*) ' WAVELENGTH = ' LAMBDA,’, FREQUENCY(FO) = ',FO 
ELSE 
WRITE(6,*) ‘ALL BOUNDARY VALUE GENERATION METHODS - DISABLED' 
DO 5, K = 1, RES 
READ(3,*) BCOND(K) 
CONTINUE 
ENDIF 


POLAR COORDINATE INPUT ROUTINE 


IF(( CHAR. EQ. 'P'). OR. (CHAR. EQ. 'p')) THEN 
READ(3,1020) THETA, RADIUS 
SFAC = 2*PI*DD/RADIUS 
MESH(0,4) = (RADIUS*SIN( THETA*P1/180.0))*SFAC + XORIGIN 
MESH(0,5) = (RADIUS*COS( THETA*PI/180.0))*SFAC + YORIGIN 
DO 10, J = 1, RES-1 
READ(3,1020) THETA, RADIUS 
MESH(J,4) = (RADIUS*SIN( THETA*P1/180.0))**SFAC + XORIGIN 
MESH(J,5) = (RADIUS*COS( THETA*P1/180.0))**SFAC + YORIGIN 
CONTINUE 
WRITE(6,*) ‘POLAR COORDINATE INPUT SELECTED ' 


RECTANGULAR COORDINATE INPUT ROUTINE 


ELSEIF( (CHAR. EQ. 'R'). OR. (CHAR. EQ. 'r')) THEN 
READ(3,1010) MESH(0,4), MESH(0,5) 
SFAC = 2*PI*DD/((MESH(0,4)**2+MESH(0,5)**2)**0. 5) 
MESH(0,4) = MESH(0,4)*SFAC + XORIGIN 
MESH(0,5) = MESH(0,5)*SFAC + YORIGIN 
hOe2o. J = 1, RES=-1 
READ(3,1010) MESH(J,4), MESH(J,5) 


MESH(J,4) = MESH(J,4)*SFAC + XORIGIN 
MESH(J,5) = MESH(J,5)*SFAC + YORIGIN 
CONTINUE 


WRITE(6,*) ‘RECTANGULAR COORDINATE INPUT SELECTED ' 
ELSE 
WRITE(6,*) ‘INPUT DATA FILE COORDINATE SPECIFICATION ERROR’ 
ENDIF 
WRITE(*,1100) DD, SFAC 
MESH(RES ,4) = MESH(0,4) 
MESH(RES,5) = MESH(0,5) 


NPNTS = J 
WRITE(6,*) "NODE AT 0 DEGREES (X/Y COORDINATES) = ', MESH(0,4), 
CMESH(0,5) 

WRITE(6,*) 'X, Y OFFSETS = ', XORIGIN, YORIGIN 


WRITE(6,1090) ALPHA, BETA 
MRES = MRESW**2*PI 

WRITE(6,**) 'MESH RESOLUTION 

DIST = DISTW*2**PI 

WRITE(6,%*) ‘CONTOUR DISTANCE 
WRITE(6,**) ‘NUMBER OF DATA POINTS 


" \MRESW,' WAVELENGTHS ' 


" ,DISTW,' WAVELENGTHS ' 
" ,RES 
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C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


C 


IF( (METHOD. EQ. 3). OR. (METHOD. EQ. 6)) THEN 
WRITE(6,*) 'START NODE 
WRITE(6,*) ‘STOP NODE 

ENDIF 


WRITE(4,*) '17' 


WRITE(4,*) ‘10 if 
55 le 
WRITE(4,*) '2 1 
80 he 
WRITE(4,*) ‘2 1 
90 Oe 
WRITE @Apc Oe a2 il 
100 ie 
WRITE(4,*) '2 1 
1G yi 
WRITE( Gcenm 2 1 
120 iL 
WRITE(4,*) '2 il 
130 Ds 
WRITE(4,*) '2 1 
140 i 
WRITE(4,*) ‘2 il 
150 ae 
WRITE(4,**) ‘2 1 
160 i 
WRITE(4,*) ‘2 a 
170 Wh! 
WRITE(4,*) '2 1 
180 1' 
WRITE(4,*) ‘2 1 
190 28 
WRT te Gece i 
200 i 
WRITE(4,*) '2 1 
210 2 
WRITE(4,*) '2 1 
220 ee 
WRITE(4,*) ‘2 1 
220 ys 


WRITE(4,1070) NAME 
Rimes) 1 

WRITE(4,*) ‘Mesh Resolution’ 
WRITE(4,*) ‘L' 

WRITE(4,1030) MRESW 
WRITE(4,*) '‘L' 

WRITE(4,*) ‘Contour Distance’ 
WRITEC2, =) - Li 

WRITE(4,1030) DISTW 
WRITE(4,*) ‘L' 

WRITE(4,*) ‘Number of Points’ 
WRITHG G0) an le 

WRITE(4,1040) RES 

WRITECS mL: 

WRITE(4,*) 'Method' 

WRIDE( 4.°3)pe al 


' | STARTND 
" | STOPND 


138 
365 
S955 
505 
365 
365 
365 
36S 
365 
365 
365 
365 
30> 
S15 
365 
305 


365 


aqIaanqa 
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WRITE(4,1040) Method 
WRITE(4,*) 'L' 
WRITE(4 ,*) Perrin 
WRITE(4,%*) ‘L' 
WRITE(4,1050) XORIGIN 
WRITE(4,*) ky 
WRITE(4,*) 'Y Origin' 
WRITE(4,*) 'L' 
WRITE(4,1050) YORIGIN 
WRITE(4,*) ‘L' 
WRITE(4,*) ‘Alpha’ 
WRITE(4,*) 'L' 
WRITE(4,1060) ALPHA 
WRITE(4,*) ‘L' 
WRITE(4,*) 'Beta' 
WRITE(4,*) ‘L' 
WRITE(4,1060) BETA 
WRITE(4,*) ‘L' 
WRITE(4,*) ‘9! 


WersO, Il = 0, NPNIS-1 

WRITE C23. IeMESH( 11,45, MESH(I1, 5) 
CONTINUE 
WRITE( 2 ,**) MESH(O , 4), MESH(O, 5) 
WRITE(2,*) ‘999992, 999991 
WRITE(2,*) ‘999990, 999990 ' 


FORMAT( Al) 

FORMAT(10X,F12. 8,4X,F12. 8) 

MORMAT( 3X,13,5X,F12.8) 

FORMAT(F8. 6) 

FORMAT(I4) 

FORMAT( 1X,F8. 6) 

FORMAT(1X,F8.6,' - J',F8.6) 

FORMAT(A12) 

FORMAT(/) 

FORMAT( 1X, ‘ALPHA bo eee hone BETA = FS. 42x, + J 
,F6. 4) 

FORMAT( 1X,‘ DESIRED DISTANCE = ',F8.5,', SClALeeracmore—— FS. 5) 

RETURN 

END 


SUBROUTINE ROTATE(CNPNTS , METHOD , STARTND , STOPND ) 


THIS SUBROUTINE ROTATES THE SURFACE POINTS TO ALLOW FOR 
A REARRANGING OF THE START NODE (FIRST DATA POINT). 


REAL MESH(0: 200,5) 
ENTEGER 1, NPNIS, METHOD, STARIND, STOPND 
COMMON/BLK1/MESH 


IF( METHOD. EQ.6) THEN 
DOSIOy i = 0, NENIS-1 
MBSHCI ,1) = MESH( 1,4) 


30 


mao 


Ci CoC: © 36G) 
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MESH(I,2) = MESH(I,5) 
CONTINUE 


DO 20, I = STARTND-1, NPNTS-1 
MESH( 1-(STARIND-1) 54)" = MESHCE 
MESH( 1 -CSTARIND-1),5) >= MEsHGie® 
CONTINUE 


DO 30, 1 = 1.7 StAkins 
MESHC I+NPNTS -STARTND , 4) 
MESHCLENENTS— sian) 
CONTINUE 
ENDIF 
STOPND = STOPND - STARTND + 1 
IF(STOPND. GI. NPNTS) THEN 
STOPND = STOPND - NPNTS 
ELSEIECSLOPND. Ei) O thin 
STOPND = NPNTS + STOPND +1 
Blob TECslORND. EO, 0) mini 


MESH(I-1,1) 
MESH(I-1,2) 


STOPND = 1 
ELSE 
STOPND = STOPND 
ENDIF 
RETURN 
END 


SUBROUTINE BOUND(MRES ,NPNTS , METHOD, STOPND) 


THIS SUBROUTINE REORDERS THE THE SURFACE POINTS 
BASED ON THE DESIKED INPUT MESH KESOLU 1) Sei 
OBJECT IS ALSO BISECTED. 


REAL PERIM, DIST, MESH(0:200,5), TEMP, TEMP1, MRES, MRESN, MRESLW 
REAL MRESNL, MRESNR, TEMPR, TEMPL, DISTB, DZ, PI, MRESW, MRESRW 
INTEGER I, PERND, BIND, NPNTS,STARTND, STOPND, METHOD, NEWND, J 
INTEGER M, L 

COMMON /BLK1/MESH 

COMMON /BLK2/PERND, BIND 


PI = 4. O*ATAN(1. 0) 
PERIM = 0.0 


DO 10, I = 0, NPNTS-1 
DIST = SQRT( (MESH(I+1,4)-MESH(I,4))**2+(MESH(I+1,5) - 
MESH(I,5) )****2) 
PERIM = PERIM + DIST 
MESH(I+1,3) = PERIM 
CONTINUE 
WRITE(6,**) ‘PERIMETER LENGTH = ', PERIM 
PERND = NINT(PERIM/MRES - AMOD( PERIM/MRES, 2. 0)) 
BIND = (PERND - 2)/2 
WRITE(6,*) 'PERIMETER NODE # = ', PERND 
MRESW = MRES/(2*PI) 
WRITE(6,%*) ‘YOUR REQUESTED MESH RESOLUTION OF ',MRESW 
IF( (METHOD. EQ. 3). OR. (METHOD. EQ. 4). OR. (METHOD. EQ. 5). OR. 
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CCMETHOD. EQ. 6)) THEN 


MRESNL 
MRESLW 


(PERIM -MESH(STOPND-1,3))/FLOAT( PERND/2) 
MRESNL/(2**PI) 

MRESNR = MESH(STOPND-1,3)/FLOAT( PERND/2) 

MRESRW = MRESNR/(2*PI) 

WRITE(6,*) '. . . HAS BEEN MODIFIED TO... ' 

WRITE(6,*) 'LEFT SIDE OF SEGMENT... ',MRESLW 

WRITE(6,*) 'RIGHT SIDE OF SEGMENT . . . ',MRESRW 
ELSE 

MRESN = PERIM/FLOAT(PERND) 

MRESW = MRESN/(2*PI) 

WRITE(6,*) '. . . HAS BEEN MODIFIED TO... *',MRESW 
ENDIF 


PERIMETER NODE INITIALIZATION (X,Y COORD) 


IF( (METHOD. EQ. 1). OR. (METHOD. EQ. 2)) THEN 
DO 30, I = 0, PERND-1 
TEMP = MRESN*I 
J=1 
IF(TEMP. GT. MESH(J,3)) THEN 
T= 741 
GOTO 20 
ENDIF 
TEMP1 = (TEMP - MESH(J-1,3))/(SQRT((MESH(J,4)-MESH(J-1,4)) 
c 2 + (MESH(J,5) - MESH(J-1,5))**2)) 
MESH(I+1,1) = MESH(J-1,4) + TEMP1**(MESH(J,4)-MESH(J-1,4)) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1*(MESH(J,5)-MESH(J-1,5)) 
CONTINUE 


ELSE 
HOso0. 1)— 0, PERND=1 
TEMPR = MRESNR*I 
TENE L SePERIM = MRESNIE(PERND - 1) 
Pa) EMER LE. MESHCSIOEND= les) THEN 


J=1 
IF(TEMPR. GT. MESH(J,3)) THEN 
J=J+1 
GOTO 40 
ENDIF 
TEMP1 = (TEMPR - MESH(J-1,3))/(SQRT((MESH(J,4)-MESH(J-1,4)) 
C %2 + (MESH(J,5) - MESH(J-1,5))%**2)) 


MESH(I+1,1) = MESH(J-1,4) + TEMP1*(MESH(J,4)-MESH(J-1,4)) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1**(MESH(J,5)-MESH(J-1,5)) 
ELSE 
J = 1 
IF(TEMPL. GT. MESH(J,3)) THEN 
Jy Sy i 
GOTO 50 
ENDIF 
TEMP1 = (TEMPL - MESH(J-1,3))/(SQRT((MESH(J,4)-MESH(J-1,4)) 
é wee2 + (MESH(J,5) - MESH(J-1,5))**2)) 
MESH(I+1,1) = MESH(J-1,4) + TEMP1*(MESH(J,4)-MESH(J-1,4)) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1*(MESH(J,5)-MESH(J-1,5)) 


CGI Ge) 


90 


100 


110 


ENDIF 
CONTINUE 
ENDIF 


BISECTION NODE INITIALIZATION (X,Y COORD) 


WRITE(6,*) 'BISECTION NODE # = ' BIND 
IF( (METHOD. EQ. 1). OR. (METHOD. EQ. 3)) THEN 
DO 70, I = 1, BIND 
TEMP1 = I/FLOAT(BIND + 1) 
MESH(I+PERND,1) = MESH(1,1) + TEMP1*(MESH(BIND+2,1) - 
CMESH(1,1)) 
MESH(I+PERND,2) = MESH(1,2) + TEMP1*(MESH(BIND+2,2) - 
CMESH(1,2)) 
CONTINUE 


ELSE 
DO 80, I = 2, PERND+1 
MESH(PERND+I-1,1) = MESH(PERND+2-I,1) + 0.5*(MESH(I,1) - 


C MESHOCPERNDE2-151)) 
MESH( PERND+1-1,2) = MESHCPERND+2-1,2) + 0. 5*( MESH Gis 
C MESH PERND=2-1525) 
CONTINUE 
ENDIF 


IF( (METHOD. EQ. 5). OR. (METHOD. EQ.6)) THEN 
DO 90, M = 1, BIND 
MESH(M,4) = MESH(PERND+M, 1) 
MESH(M,5) = MESH(PERND+M, 2) 
CONTINUE 
Dicive— ono 
MESH(0,3) = 
DO 100, L= BIND+1 
IF(L. EQ. 1) THEN 
DIST = SQRT((MESH(1,1) - MESH(1,4))**2 + 
C(MESH(1,2) - MESH(1,5))**2) 
ELSEIF(L. EQ. BIND+1) THEN 
DIST = SQRT((MESH(BIND,4) - MESH(BIND+2,1)) 
C*"2 + (MESH(BIND,5) - MESH(BIND+2,2) )**2) 
ELSE 
DIST = SQRT((MESH(L-1,4) - MESH(L,4)) 
C**2 + (MESH(L-1,5) - MESH(L,5))**2) 
ENDIF 
DISTB = DISTB + DIST 
MESH(L,3) = DISTB 
CONTINUE 
DZ = DISTB/(BIND + 1.0) 
WRITE(6,*) 'DZ SPACING = ee 


O20 
1, 


DO 120, I 
TEMP 
J=1 
IF(TEMP. GT. MESH(J,3)) THEN 


Pb LIND 
D2*I 
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yy 
GOTO 110 
ENDIF 
IF(J.EQ.1) THEN 
TEMP1 = TEMP/(SQRT((MESH(1,1) - MESH(1,4))**2 + 
C (MESH(1,2) - MESH(1,5))**2)) 
MESH( PERND+1,1) = MESH(1,1) + TEMP1*(MESH(1,4) 
é - MESH(1,1)) 
MESH( PERND+1,2) = MESH(1,2) + TEMP1*(MESH(1,5) 
C - MESH(1,2)) 
ELSE 
TEMP1 = (TEMP - MESH(J-1,3))/(SQRT((MESH(J,4) - 
C MESH(J-1,4))**2 + (MESH(J,5) - MESH(J-1,5))**2)) 
MESH( PERND+1,1) = MESH(J-1,4) + TEMP1*(MESH(J,4) 
C - MESH(J-1,4)) 
MESH( PERND+1,2) = MESH(J-1,5) + TEMP1*(MESH(J,5) 
c - MESH(J-1,5)) 
ENDIF 
120 CONTINUE 
ENDIF 
RETURN 
END 


i a @ 


SUBROUTINE NORMAL 


THIS ROUTINE COMPUTES THE X AND Y COMPONENTS 
OF THE OUTWARD UNIT NORMAL AT EACH SURFACE POINT. 


SP al SP Tae ps aD 


REAL DR, DZ, DL, MESH(0:200,5) 
INTEGER I, PERND, BIND 
COMMON/BLK1/MESH 
COMMON/BLK2/PERND, BIND 
DO 10, I = 1, PERND 
IF(I.EQ.1) THEN 
DR = MESH(2,1) - MESH(PERND,1) 
DZ = MESH(2,2) - MESH(PERND,2) 
DL = SQRT(DR*DR+DZ*""DZ) 
MESH(1,3)=-DZ/DL 
MESH(1,4)=DR/DL 
ELSEIF(I.EQ. PERND) THEN 
DR = MESH(1,1) - MESH(PERND-1,1) 
DZ = MESH(1,2) - MESH(PERND-1,2) 
DL = SQRT(DR*“DR+DZ*DZ) 
MESH( PERND ,3)=-DZ/DL 
MESH( PERND ,4)=DR/DL 
ELSE 


MESH(I+1,1) - MESH(I-1, 1) 
MESH(I+1,2)- MESH(I-1,2) 
SQRT( DR**DR+DZ*DZ) 
MESH(1 ,3)=-DZ/DL 
MESH(1,4)=DR/DL 
ENDIF 
10 CONTINUE 
RETURN 


ww) 
N 
| 
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END 


SUBROUTINE NODSET( METHOD , NABC , NBMX , UNK) 


THIS ROUTINE USES A TWO-SWEEP TECHNEGQUES IO COM Ul thn 
NUMBER OF NODES ALONG EACH NODE ROW, I. A MESH 
ORIENTATION ATTRIBUTE IS ALSO SET. 


SET Z-AXIS SPACING AND ENDPOINT NODES 


REAL DZ, 2Z, MESH(0:200,5), D, DIST, DISTB 

INTEGER I, J, NODES(200), MOR(200), PERND, NOLD, NNEW, BIND 
INTEGER L, METHOD, NABC(100,3), NBMX, UNK 

CHARACTER*1 CHAR2, CHAR3, CHAR4 

COMMON/BLK1/MESH 

COMMON /BLK2/PERND , BIND 

COMMON/BLK3/NODES 

COMMON /BLK6/MOR 

COMMON/BLK9 /CHAR2 , CHAR3 , CHARS 


UNKe—a0 
NBMX = -99 
IF( (METHOD. EQ. 1). OR. (METHOD. EQ. 3)) THEN 
DZ = (SQRT((MESH(1,2) - MESH(BIND+2,2))%**2 + 
e (MESH(1,1) - MESH(BIND+2,1))***2))/(BIND + 1.0) 
ELSE 
DISTB = 0.0 
DO 10, L = 1, BIND#+1 
IF(L. EQ. 1) THEN 
DIST = SQRT((MESH(1,1) - MESH(PERND+1,1))**2 + 
C(MESH(1,2) - MESH( PERND+1, 2) )****2) 
ELSEIF(L. EQ. BIND+1) THEN 
DIST = SQRT((MESH( PERND+BIND,1) - MESH(BIND+2,1)) 
C**2 + (MESH(PERND+BIND,2) - MESH( BIND+2,2))***2) 
ELSE 
DIST = SQRT((MESH( PERND+I-1,1) - MESH( PERND+I,1)) 
C2 + (MESH(PERND+I-1,2) - MESH( PERND+I ,2))**2) 
ENDIF 
DISTB = DISTB + DIST 
CONTINUE 
DZ = DISTB/(BIND + 1.0) 
ENDIF 
WRITE(6,*) 'BISECTION SPACING =', DZ 
NODES(1) = 2 
NODES(2) = 3 
NODES(BIND+2) = 2 
NODES( PERND+1) = 2 
NODES( PERND) = 3 
PERFORMING FORWARD-SWEEP 
DO 20 I = 3, BIND+1 
NOLD = NODES(I-1) 
D = SQRT((MESH(I,1) - MESH(I+PERND-1,1))**2 + 
C (MESH(1,2)-MESH( I+PERND-1,2))**2) 
NNEW = INT(0.1 + D/DZ) + 2 
NODES(I) = NNEW 
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IF (NNEW. GT. NOLD+1) NODES(I) = NOLD + 

IF (NNEW. LT. NOLD-1) NODES(I) = NOLD - 

Pee CNOUNS Cl). Uline NODESE)) =.35 
CONTINUE 


> 


J = PERND + 2 
DO 30, I = PERND-1, BIND+3, -1 
NOLD = NODES(I+1) 
D = SQRT((MESH(I,1) - MESH(J,1))%#**2 + 


C CUESHCIs 2) -MESH J+ 2ga22)) 


NNEW = INT(0.1 + D/DZ) + 2 
NODES(1I) = NNEW 
IF (NNEW. GT. NOLD+1) NODES(I) 
IF (NNEW. LT. NOLD-1) NODES(I) 
IF (NODES(I).LT.3) NODES(I) = 3 
J=J+1 

CONTINUE 


NOLD + 
NOED = 1 


— 


BACK-SWEEP TO RESET LAST NODES IF NEEDED 


I BND 2 

Ai Tey 

Pee iwe.2) GO TO 50 

Pee NOBESC 1). LE. NODES(I+1)+1) GO TO 60 

ewes (i) = NODES(I+1) + 1 

GO TO 40 

CONTINUE 

WRITE(6,*) ' PROGRAM ABORTED IN NODSET RIGHT SIDE BACKSWEEP ' 
SOP 


CONTINUE 
im oIND + 2 
m= 1+ 1 


feeCi. EO.PERND) GO TO 80 

MONOURS CT). LE. NODES(I-1)+1) GO TG 90 

NObHS( i) = NODES(I-1) + 1 

GO TO 70 

CONTINUE 

WRITE(6,*) ' PROGRAM ABORTED IN NODSET LEFT HALF BACKSWEEP ' 
STOP 

CONTINUE 


FORWARD SWEEP TO LOAD MESH ORIENTATION ARRAY, MOR 


MOR(1) = 0 
MOR(BIND+2) = 0 
MOR( PERND+1) = 0 


On TOO ee i=2, BIND+] 

IF(NODES( 1I+1).GT.NODES(1)) THEN 
MOR(I) = 0 

ELSEIF(NODES(I+1). LT.NODES(I)) THEN 
MOR(I) = 1 

ELSE 
IF(NODES( I+2).GT.NODES(I)) THEN 

MOR(I) = 0 
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ELSEIF(NODES(1I+2). LT. NODES(I)) THEN 
MOR(I) = 1 
ELSE 
MOR(I) = MOR(I-1) 
ENDIF 
ENDIF 
100 CONTINUE 


DO 110, I = PERND, BIND+3, -1 
IF(NODES(I-1).GT. NODES(I)) THEN 
MOR(I) = 0 
ELSEIF(NODES(I-1). LT. NODES(1I)) THEN 
MOR(I) = 1 
ELSE 
IF(NODES(I-2).GT. NODES(1I)) THEN 
MOR(I) = 0 
ELSEIF(NODES(I-2). LT. NODES(I)) THEN 
MOR(I) = 1 


MOR(I) 
ENDIF 
ENDIF 
110 CONTINUE 
C LOAD NABC ARRAY 
DO 120, I = 1, BIND+2 
IF(I.EQ.1) THEN 
NABC(1,1) = 
NABC(1,2) = 


MOR( I+1) 


Wr © 


ELSEIF(I. LE. BIND) THEN 
NABC(I,1) = NABC(I-1,2) 
NABC(1,2) = NABC(I-1,3) 
NABC(I,3) = NODES(PERND-I+1) + NODES(I+1) - 3 
ELSEIF(I.EQ. BIND+1) THEN 
NABC(I,1) = NABC(I-1,2) 
NABC(I,2) = NABC(I-1,3) 
NABC(I,3) = 1 
ELSE 
NABC(I,1) 
NABC(I,2) 
NABC(I,3) 
ENDIF 
120 CONTINUE 


NABC( I-1, 2) 
NABC(I-1,3) 
0 


DO 130, I = 1, BIND+2 

IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(12,**) I,NABC(I,1),NABC(I,2),NABC(I,3) 

ENDIF 

UNK = UNK + NABC(I,2) 

IF(NABC(1I,2).GE.NBMX) THEN 
NBMX = NABC(I,2) 

ENDIF 

130 CONTINUE 


WRITE(*,*) ‘MAXIMUM UNKNOWN WIDTH 
WRITE(*,*) ‘TOTAL # OF UNKNOWNS 


" /NBMX 
' , UNK 


90 
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Ca-G2) CY C2 C2 C2 


WRITE(*,*) ' 
RETURN 
END 


SUBROUTINE SORTER(1,LEL, LAND) 


THIS SUBROUTINE GENERATES A MESH ROW FOR THE INPUT ROW I. 
LOADING OF THE LOCAL NODE-ELEMENT CONNECTION MATRICES LND AND 
NDL FOR ELEMENTS BETWEEN I AND I+1 NODE ROWS. REFERENCE IS TO 
THE LEFT SIDE OF THE Ith ROW OR VECTOR. 


INTEGER NODES(200), MOR(200), LND(0:200,3), NDL(200,4), NCT(200) 

INTEGER I, PERND, BIND, LEL, LAND, J, K, LL, JJ, NN, KK, N, Nl, N2 

INTEGER NDMX, NDS1L, NDS2L, NDS1R, NDS2R, LMX 

COMMON /BLK2/PERND, BIND 

COMMON /BLK3/NODES 

COMMON /BLK4/LND ,NDL,NCT 

COMMON /BLK6/MOR 

IF(I.EQ. BIND+2) THEN 
WRITE(6,*) 'ERRORED OUT IN SUBROUTINE SORTER, YOU ATTEMPTED’ 
WRITE(6,*) ' TO CALL SORTER WITH I = BIND + 2!' 
WRITE(6,*) ' THIS ROW HAS NO ELEMENTS' 
RETURN 

ENDIF 

hero, J = 0, 200 
POO. kK = 1, 

LND(J,K) 

CONTINUE 

CONTINUE 

NDMX = 200 

NDS1L = NODES( PERND+2-1) 

NDS2L = NODES(PERND+1-I) 

NDS1R = NODES(I) 

NDS2R = NODES(I+1) 

LMX = NDSIL + NDS1R + NDS2L + NDS2R - 2 


S 
= 0 


LEFTSIDE OF THE BISECTION SEGMENT 
TOP ROW 


IF(I.EQ.1) THEN 
LND(1,1) = 
LND( 1,2) 
LND( 1,3) 
LND( 2, 1) 
LND( 2,2) 
LND( 2,3) 
LND( 3,1) 
LND( 3,2) 
LND( 3,3) 
LND( 4,1) 
LND( 4,2) 
LND( 4,3) 


bund ninbud due te wet 
ANNOOHLNYM OD Fe WOH 
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C2 C2. C2 


LND(5,1) 
LND(5 ,2) 
LND(5 ,3) 
LND( 6,1) 
LND( 6,2) 
LND(6,3) 
LEL = 6 
LAND = 7 


BOTTOM ROW 


IND(1,1) = 
LND( 1, 2) 
LND(1,3) 
LND( 2,1) 
LND( 2,2) 
LND( 2,3) 
LND(3, 1) 
LND(3, 2) 
LND(3,3) 
LND(4,1) 
LND( 4,2) 
LND( 4,3) 
LND(5,1) 
END ES) 
LND(5,3) 
LND(6,1) 
LND( 6,2) 
LND( 6,3) 


“SIDR OVD OF 


UN NEW NSW DM SID SID DAE OPO 


LEL = 


LAND = 


EQUAL NODE 


6 
7 


NUMBERS 


ELSELE CI. bO sD ee aiin 


ELSEIFC(NDS1L. EQ. NDS2L) THEN 
O (LH ORIENTATION) 
TPC MOR UPERKND-2- 4 EQ. 0) GHEN 


FOR MOR = 


END@ 11) = 
LND(1,2) = ion + NDS1R 
LND(1,3) = 2 
LND(2,1) = NDS1L + NDS1R + 1 
LND(2,2) = 2 
LND(2,3) = NDSIL + NDS1R 
DO 40, N = 1, NDS1L-2 
N1 = 2*N + 1 
N2 = Ni + 1 
DO 30, K = 1, 3 
LND(N1,K) = LND(1,K 
LND(N2,K) = LND(2,K 
CONTINUE 
CONTINUE 
1 (RH ORIENTATION) 
LND(1,1) = NDS1L + NDS1R 
MDC, 2D) =i 


LND(1, 3) 
LND( 2,1) 
LND( 2,2) 


NOS NDSIR es 1 
2 
NUS ie NDS AR 1 


50, K=1, 3 
LND(N1,K) = LND(1,K) +N 
LND(N2,K) = LND(2,K) +N 
50 CONTINUE 
60 CONTINUE 
ENDIF 


LEFTHAND MESH ORIENTATION 


Carga €> 


Pioet’( MOR( PERND+2-1).EQ.0) THEN 
LND( 1,1) NES Iie NDS ee 
NDC 1, 2) 1 
LND( 1,3) NDS1L + NDSIR 

0 

NDSIL + NDSIR 


ie 
Ff 
7) 
(aeN 
© 
_ 
ww 
ri ue 


DO 70, K = 1, 3 
LND(N1,K) = LND(O,K) + N 
LND(N2,K) = LND(1,K) +N 
70 CONTINUE 
80 CONTINUE 


C RIGHTHAND MESH ORIENTATION 


ELSE 
LND(1,1) 
LND( 1,2) 
LND(1,3) 
LND(0,1) 
LND(0,2) = 1 
LND(0,3) = NDS1L + NDS1R 
DO 100, N= 1, NDSIL-2 


2 

NDS1L + NDSI1R 

1 

Nadie NDS IR sal 


LND(N1,K) = LND(O,K 
LND(N2,K) = LND(1,K 
90 CONTINUE 
100 CONTINUE 
ENDIF 


eee 
oN 


eel Oa) ]OR.( 1. EQSBIND+1)) THEN 
COTO. 230 
ENDIF 
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Ca Git CG Care) 


>) 


LED =shz 


RIGHTSIOE OF THE BISEC TION SEGIEnNT 


EQUAL NODE NUMBERS 


IFCNDSIR. EQ. NDS2K)) THEN 
MOR = O (LH ORIENTATION) 
TF(MORCI) 2 2e7 0) Shen 


LNDGEE Eis) NDS1IL + NDOSIRS ENS 2 aa 


LND(LEL+1,2) = NDS1L 
LND(LEL+1,3) = NDS1L + NDS1R + NDS2L 
LND( LEL+2,1) = NDS1L + 1 
LND(LEL+2,2) = NDS1L + NDS1R + NDS2L 
LND(LEL+2,3) = NDS1L 
DO 120, N = 1, NDS1R-2 
Ni = 2*N + 1 + LEL 
N2 = Nl +1 
DO 110, K = 1, 3 
LND(N1,K) = LND(LEL+1,K) + N 
LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 
LEL = N2 
LAND = LMX 
MOR = 1 (RH ORIENTATION) 
ELSE 
LND(LEL+1,1) = NDSIL 
LND(LEL+1,2) = NDS1L + NDS1R + NDS2L - 1 
LND(LEL+1,3) = NDS1L + 1 
LND(LEL+2,1) = NDSIL + NDS1R + NDS2L 
LND(LEL+2,2) = NDS1L + 1 
LND(LEL+2,3) = NDS1L + NDS1R + NDS2L - 1 
DO 140, N = 1, NDS1R-2 
hie hee Loe 
N2 =N1 +1 
DO 130, K = 1, 3 
LND(N1,K) = LND(LEL+1,K) + N 
LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 
ENDIF 
LEL = N2 
LAND = LMX 


LEFT HAND MESH ORIENTATION 


ELSEIF(MOR(1I).EQ.0) THEN 
LND(LEL+1,1) = NDS1L + NDS1R + NDS2L - 1 


LND( LEL+1,2) = NDS1L 
LND( LEL+1,3) = NDS1L + NDS1R + NDS2L 
LND(LEL+2,1) = NDSIL + 1 


END Cin 2) NDSIL + °NDSIR = NDSZL 
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LND( LEL+2,3) = NDS1IL 


NDS1R + NDS2L 
NDSIR + NDSZE 


NESIRet NDS2ZL 


é 
DO 160, N = 1, NDS1R-1 
N1 = 2*N + LEL + 1 
me 180, 6 = 1, 8 
LND(N1,K) = LND(LEL+1,K) +N 
150 CONTINUE 
160 CONTINUE 
c 
DO 180, N = 1, NDS1R-2 
N2 = 2*N + LEL + 2 
DO 170, K = 1, 3 
LND(N2,K) = LND(LEL+2,K) + N 
170 CONTINUE 
180 CONTINUE 
C 
LEL = Nl 
LAND = LMX 
C 
c RIGHT HAND MESH ORIENTATION 
E 
ELSE 
LND(LEL+1,1) = NDS1L 
LND( LEL+1,2) = NDS1L + 
LND(LEL+1,3) = NDS1L + 1 
LND( LEL+2,1) = NDS1L + 
LND( LEL+2,2) = NDS1L + 1 
LND( LEL+2,3) = NDS1L + 
ic 
DO 200, N = 1, NDS1R-2 
N1 = 2*N + LEL + 1 
no U9, = hy s 
LND(N1,K) = LND(LEL+1,K) + 
190 CONTINUE 
200 CONTINUE 
C 
DO 220, N = 1, NDS1R-3 
Ne = 27N 2° LEn + 2 
DO 210, K = 1, 3 
LND(N2,K) = LND(LEL+2,K) + 
210 CONTINUE 
220 CONTINUE 
C 
LEL = Nl 
LAND = LMX 
C 
C 
ENDIF 
C 
C 
C 
c ZEROING NDL AND NCT PRIOR TO FILL 
> 
230 DO 250, N= 1, NDMX 
NCT(N) = 0 
DO 240, K=1, 4 
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NDL(N,K) = 0 
CONTINUE 
CONTINUE 


SCANNING LND ARRAY TO LOAD NDL 


DO 270 LL = 1, LEL 
DO 260 JJ =1, 3 
NN = LND(LL,JJ) 
NCT(NN) = NCT(NN) + 1 
KK = NCT(NN) 
NDL(NN,KK) = LL 
CONTINUE 
CONTINUE 


END 


SUBROUTINE FINDER(I, LAND ,DIST) 


THIS SUBROUTINE DETERMINES THE (X,Y) COORDINATES OF EACH 
NODE IN THE CALLING Ith ROW OR VECTOR MESH. 


INTEGER I, J, LAND, PERND, BIND, NODES(200) 
REAL MESH(0: 200,5), NDP(200,2), DIST 
COMMON /BLK1/MESH 
COMMON /BLK2/PERND , BIND 
COMMON/BLK3/NODES 
COMMON /BLKS /NDP 
IF(I.EQ.1) THEN 
NDP(1,1) = MESH(1,1)+MESH(1,3)*DIST 


NDP(1,2) = MESH(1,2)+MESH(1,4)*DIST 
NDP(2,1) = MESH(1,1) 

NDP(2,2) = MESH(1,2) 

NDP(3,1) = MESH(PERND, 1)+MESH(PERND ,3)*DIST 
NDP(3,2) = MESH(PERND, 2)+MESH(PERND ,4)**DIST 
NDP(4,1) = MESH(PERND,1) 

NDP(4,2) = MESH(PERND,2) 

NDP(5,1) = MESH( PERND+1,1) 

NDP(5,2) = MESH(PERND+1,2) 

NDP(6,1) = MESH(2,1) 

NDP(6,2) = MESH(2,2) 

NDP(7,1) = MESH(2,1)+MESH(2,3)*DIST 
NDP(7,2) = MESH(2,2)+MESH(2,4)*DIST 


ELSEIFC1. EQ. BIND+1) THEM 
NDP( 1,1) = MESH(BIND+3, 1)+MESH( BIND+3 ,3)*DIST 


NDP(1,2) = MESH( BIND+3 ,2)+MESH(BIND+3,4)*DIST 
NDP(2,1) = MESH(BIND+3, 1) 

NDP(2,2) = MESH( BIND+3, 2) 

NDP(3,1) = MESH(PERND+BIND, 1) 

NDP(3,2) = MESH( PERND+BIND,2) 

NDP(4,1) = MESH(BIND+1, 1) 

NDP(4,2) = MESH( BIND+1, 2) 

NDP(5,1) = MESH( BIND+1,1)+MESH(BIND+1,3)*DIST 
NDP(5,2) = MESH( BIND+1,2)+MESH(BIND+1 ,4)*DIST 
NDP(6,1) = MESH( BIND+2, 1)+MESH(BIND+2,3)*DIST 
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RIESE 1 


ELSE 


ENDIF 
RETUR 


NDP(6,2) = MESH(BIND+2 ,2)+MESH( BIND+2,4)*DIST 
NDP(7,1) = MESH(BIND+2,1) 

NDP(7,2) = MESH(BIND+2,2) 

F(I. EQ. BIND+2) THEN 


WRITE(6,*) ' ERRORED OUT IN SUBROUTINE FINDER, YOU ATTEMPTED' 


tou il 


WRITE(6,*)  ' TO CALL FINDER WITH I = BIND + 2!' 
WRITE(6,*) ' THIS ROW HAS NO ELEMENTS AND NO COORDINATES' 
NODE 1 

NDP( 1,1) = MESH( PERND-I+2,1)+MESH( PERND-I+2 ,3)*DIST 


NDP(1,2) = MESH( PERND-I+2,2)+MESH( PERND-1I+2,4)*DIST 

NODE 2 TO THE BISECTION SEGMENT 

DO 10, J = 2, NODES(PERND-I+2) 
NDP( J, 1)=MESH( PERND-I+2,1)-(J-2)**(MESH( PERND-I+2,1)- 
MESH( PERND+I-1,1))/(NODES( PERND-I+2) -2) 
NDP( J, 2)=MESH( PERND-I+2 ,2)-(J-2)*(MESH( PERND-1I+2,2)- 
MESH( PERND+1-1, 2) )/(NODES( PERND-I+2) -2) 

CONTINUE 

BISECTION SEGMENT TO THE RIGHTSIDE SURFACE 

pen20)) J = Se NOnEs(l) 
NDP( NODES( PERND-1+2)+J-2,1)=MESH( PERND+I-1,1)+(J-2)* 
(MESH( 1, 1)-MESH( PERND+I-1,1))/(NODES(I)-2) 
NDP( NODES( PERND-1I+2)+J-2,2)=MESH( PERND+I-1,2)+(J-2)* 
(MESH( 1 ,2)-MESH( PERND+I-1,2))/(NODES(I)-2) 

CONTINUE 

Ith ROWS LAST NODE 

NDP( NODES( PERND-1+2)+NODES(1I)-1,1)=MESH(I,1)+MESH(1I,3)*DIST 

NDP( NODES( PERND-1+2)+NODES(1)-1,2)=MESH(I,2)+MESH(1I,4)*DIST 

I+lth ROWS FIRST NODE 

NDP( NODES( PERND-1+2)+NODES(1),1)=MESH( PERND-I+1,1)+MESH 

(PERND-I+1,3)*DIST 

NDP( NODES( PERND-I+2)+NODES( 1) ,2)=MESH( PERND-I+1 ,2)+MESH 

(PERND-I+1,4)*DIST 

I+1TH ROW (NODE 2) TO THE BISECTION SEGMENT 

DO 30, J = 2, NODES(PERND-I+1) 
NDP( J+NODES( PERND-1+2)+NODES(1I)-1, 1)=MESH( PERND-I+1, 
1)-( J-2)**(MESH( PERND-I+1, 1)-MESH( PERND+I,1))/ 
(NODES( PERND-I+1)-2) 
NDP( JtNODES( PERND-1I+2)+NODES(1I)-1,2)=MESH( PERND-I+1, 
2)-(J-2)*(MESH( PERND-I+1,2)-MESH( PERND+I ,2))/ 
(NODES( PERND-I+1)-2) 

CONTINUE 

I+1th ROW BISECTION SEGMENT TO THE RIGHTSIDE SURFACE 

DO 40, J = 3, NODES(I+1) 
NDP( LAND-NODES(I+1)+J-1,1)=MESH(PERND+I1 ,1)+(J-2)* 
(MESH( I+1,1)-MESH( PERND+I , 1) )/(NODES(I+1)-2) 
NDP( LAND-NODES( I+1)+J-1,2)=MESH(PERND+I1 ,2)+(J-2)* 
(MESH( I+1,2)-MESH( PERND+I , 2) )/(NODES( I+1)-2) 

CONTINUE 


LAST NODE 

NDP(LAND,1) = MESH(I+1,1)+MESH(I+1,3)*DIST 
NDP(LAND,2) = MESH(I+1,2)+MESH(1I+1,4)*DIST 
N 


OF, 


LS 1 Bp 


C2 C2 C2: C25) Gi GCG 0 6) CG Gare) 


END 


SUBROUTINE VARINT(J,F,ALPHA, BETA, AREA, LND) 


GENERATING VARIATIONAL FINITE ELEMENT AREA INTEGRATIONS OF THE 
LINEAR BASIS FUNCTION LAGRANGIAN FOR THE HELMHOLTZ EQUATION. 
THESE ARE RETURNED IN F(3,3). X(€3) AND Y(3) ARE THE WAVENUMBER 
NORMALIZED CARTESIAN COORDINATES OF THE TRIANGLE VERTICES. 

X = Ko*x, Y = Ko*y ALPHA AND BETA ARE COMPLEX 

MATERIAL PARAMETERS WITHIN THE ELEMENT. 


FOR TM INCIDENCE: ALPHA 
FOR TE INCIDENCE: ALPHA 


Tuc; SBE TA 
l/fer; BETA 


er 
ur 


OUTPUIS= F(3,3) - FINITE ELEMENT AREA INTEGRATION 
AREA - AREA OFF A TRIANGLE 


COMPLEX ALPHA, BETA, B12, F(3,3) 

REAL NDP(200,2), X(3), Y(3), T(3,3), AREA, DET 
INTEGER L, K, J, LND(0: 200,3) 

COMMON/BLK5/NDP 


X(1) NDPCGENDCI a) 


X(2) = NURGINDC 2) 
X(3) = NDP( LND(J,3),1) 
Y(1) = NDP(LND(J,1),2) 
¥(2) = NDP(LND(J,2),2) 
Y(3) = NDP(LND(J,3),2) 
DET = X(2)*Y(3) + X(3)*Y(1) + X(1)*Y(2) - X(3)*Y(2) - 


CXC 1) G3 2 a) 


AREA = ABS(0. 5**DET) 
B12 = BETA/12. 


T(1, 1) = (12 ee Be 

T(1,2) = (X03) ee yy a 

TC1,3)) = (1) Seal py ean 

T( 2,1) = (AC C2 ee 

TC 2 2) AO ee) Ee 

T0252) = X02) Pe 

TC3,, 1) = "CAC 2) sae) 2) Ee 
T0352) = CXC 3) “VCD See e02))) Par 
TC35 3) = (ACT) 02 ee ) et) be 
DO. AO a= 


Do! 10. b= 
F(K,L) = ALPHA*(T(1,K)*T(1,L) + T(2,K)*T(2,L)) - B12 
IF(K.EQ.L) F(K,L) = FCK,L) - B12 

F(K,L) = AREA*F(K,L) 


RETURN 
END 


SUBROUTINE LOADER(BCOND ,OFFSET,ALPHA ,BETA,NABC,IMX,NBMX,EO,SURBC, 
CCHAR1, LINE,MODE , XORIGIN, YORIGIN, SVEC , CHARS ) 


COMPLEX A(50,50),B(50,50) ,C(50,50) ,P(50, 100) 
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Pe 


50 


COMPLEX F(3,3),FROW(100,3,3),BCOND(100) , LINE(50) 

COMPLEX ALPHA, BETA, DETERM, SURBC(100), SVEC(50,100) 

REAL OFFSET ,MINAREA,MAXAREA, AREA, RATIO, E0,KO,XORIGIN, YORIGIN 
REAL UBCOND 

INTEGER 1,J,JD,K,L,NDTOP,NDBOT,NDTOT,NOD,KND,ND(3),LEL, LAND 
INTEGER LND(0: 200,3), NDL(200,4), NCT(200), PERND, BIND, JJ 
INTEGER NODES(200), MINROW, MINEL, MAXROW, MAXEL, TCALL, NL 
INTEGER N,M,NBMX,NABC(100,3),INORM, IMX,MODE 

CHARACTER*1 CHAR1, CHAR2, CHAR3, CHAR4, CHARS 


COMMON/BLK2/PERND , BIND 

COMMON/BLK3/NODES 

COMMON /BLK4/LND ,NDL,NCT 

COMMON /BLK7 /MINAREA , MINROW , MINEL, MAXAREA , MAXROW , MAXEL, AREA 
COMMON/BLK8/A,B,C,P 

COMMON/BLK9/CHAR2, CHAR3, CHAR4 


UBCOND = 1.0 
INORM = 0 

IMX = BIND + 2 
T= 1 


TCALL = BIND + 1 


IF((CHAR3. EQ. 'D’ ). OR. (CHAR3. EQ. 'd')) THEN 
WRITE( 20,1030) 
WRITE( 20,1040) 
WRITE( 20,1050) 
WRITE( 20,1060) 
WRITE( 20,1070) 
ENDIF 
WRITE(*,1000) I,TCALL 


CALL SORTER(I,LEL, LAND) 

CALL FINDER(I,LAND,OFFSET) 

LF(. NOT. ((CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 
CALL BNDC(1,BCOND,EO,SURBC,CHAR1 , ALPHA, BETA, LINE, MODE, XORIGIN, 


CYORIGIN, CHARS) 


ENDIF 

CALL FILL(1,LEL,FROW, ALPHA, BETA) 

IF(. NOT. ((CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 

Sc ZERO 

ENDIF 

ESTABLISH THE NUMBER OF NODES IN THE I = 1 AND I = 2 ROWS 


NDTOP = 2 
NDBOT = 5 
NDTOT = 7 


Ene PF FOR I = 1 (TOP) 


J=1 
fo70. JD = 1, NCTC2) 
L = NDL(2,JD) 
pO 50, K=1, 3 
ND(K) = LND(L,K) 
IF(ND(K).EQ. 2) KND = K 
CONTINUE 


og 


DO 60 0K, =i 
NL = ND(K) 
C BOUNDARY CONDITION 
TPCNT EOS et iia 
P(J,1) = -UBCOND*FROW(L,KND,K) + P(J,1) 


e UPPER ROW 
ELSEIF(NL. EQ. 2) THEN 
B(J,1) = FROW(L,KND,K) + B(J,1) 
@ LOWER ROW 
ELSEIF( (NL. GT. 3). AND. (NL. LT. 7)) THEN 
C(J,NL-3) = FROW(L,KND,K) + C(J,NL-3) 
C ERROR 
ELSE 
WRITE(*,**) NL,'ERROR - INDEX EQUALS 3 OR 7' 
ENDIF 
60 CONTINUE 
70 CONTINUE 
C 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(1,*) I 
DO 72. M = eee > 
WRITE(1,1020) (REAL(B(M,N)), N 
72 CONTINUE 
WRG Cl. oie ee 
DOGeelke —=2 i a NeRee 
WRITE(1,1020) (REAL(C(M,N)), N 
73 CONTINUE 
WRITE(1,*) ' ' 
DO 74.0 ie=ai NABCOM. 2) 
WRITE(1,1020) (REAL(P(M,N)), N = 1, PERND) 
74 CONTINUE 


1, NABC(I,2)) 


iH 


1, NABC(I,3)) 


¢ 
WRITE(1,*) ' ' 
WR IE (oles -o)ae 
WRITE@IG- am 
ENDIF 
C 
IF(. NOT. ((CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 
CALL MARCH(1I, IMX,NBMX,NABC(I,1),NABC(1I,2),NABC(1,3),INORM,SVEC) 
CALL. ZERO 
ENDIF 
6 BEGIN LOOPING 
DO 900, I = 1, BIND 
6 
DO 400, J = 1, NDBOT-2 
NOD = J + NDTOP + 1 
DO 300, JD = 1, NCTCNOD) 
L = NDL(NOD, JD) 
DO 100, K = 1, 3 
ND(K) = LND(L,K) 
IF(ND(K).EQ. NOD) KND = K 
100 CONTINUE 
DO 200, K=1, 3 
NL = NDCK) 
C BOUNDARY CONDITION 


[FC CL7EQ. 1) CANDSt NES EO] 1) me tneN 


100 


P(J,NL) = -UBCOND*FROW(L,KND,K) + P(J,NL) 
ELSEIF(NL. EQ. 1) THEN 
P( J, PERND-1I+2 )=-UBCOND**FROW( L,KND ,K)+P(J, PERND-I+2) 
@ SPECIAL CASE FOR NODE #2 AND I = 1 
ELSEIF((I. EQ. 1). AND. (NL. EQ. 2) )THEN 
A(J,1) = FROW(L,KND,K) + A(J,1) 
ELSEIF(NL. EQ. NDTOP) THEN 
P(J,1) = -UBCOND*FROW(L,KND,K) + P(J,I) 
ELSEIF(NL. EQ. NDTOP+1) THEN 
P(.J, PERND-I+1)=-UBCOND*FROW(L,KND,K)+P(J,PERND-I+1) 
ELSEIF(NL. EQ. NDTOT) THEN 
P(J,I+1) = -UBCOND*FROW(L,KND,K) + P(J,I+1) 
C UPPER ROW 
ELSEIF(NL. LT. NDTOP) THEN 
A(J,NL-1) = FROW(L,KND,K) + AC(J,NL-1) 
C LOWER ROW 
ELSEIF(NL. LT. NDTOT) THEN 
B(J,NL-NDTOP-1) = FROW(L,KND,K) + B(J,NL-NDTOP-1) 
t ERROR 


ELSE 
WRITE(*,**) NL,'ERROR - DUE TO INDEX GREATER THAN NODE TOTAL ' 
ENDIF 
e 
200 CONTINUE 
300 CONTINUE 
400 CONTINUE 
c 
9 INDEX FOR NEXT ROW AND COMPLETE B, C, P FILLS 
g 
JJ =I+1 
WRITE(6,1010) JJ,TCALL 
CALL SORTER(JJ,LEL, LAND) 
CALL FINDER(JJ,LAND , OFFSET) 
IF(. NOT. (( CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 
CALL BNDC(JJ,BCOND,E0O,SURBC,CHAR1,ALPHA, BETA, LINE,MODE, 
3 XORIGIN , YORIGIN ,CHAR4) 
ENDIF 
CALL FILL( JJ, LEL,FROW, ALPHA, BETA) 
C ESTABLISH THE NUMBER OF NODES IN THE I+1th AND I+2th ROWS 
IF(JJ. NE. (BIND+1)) THEN 
NDTOP = NODES( PERND+2-JJ) + NODES(JJ) - 1 
NDBOT = NODES(PERND+1-JJ) + NODES(JJ+1) - 1 
NDTOT = NDBOT + NDTOP 
ELSE 
NDTOP = 5 
NDBOT = 2 
NDTOT = 7 
ENDIF 
‘@ 


DO 800, J = 1, NDTOP-2 
NOD = J +1 
DO 700, JD = 1, NCT(NOD) 
L = NDL(NOD, JD) 
© HOO, 1K Sal, s 
ND(K) = LND(L,K) 
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IF(ND(K). EQ: NOD) KND = K 


500 CONTINUE 
DO 600, K = 1, 3 
NL = ND(K) 
E BOUNDARY CONDITION 


IF(NL. EQ. 1) THEN 
P( J, PERND-JJ+2 )=-UBCOND*FROW(L,KND ,K)+P(J, PERND-JJ+2) 
ELSEIF(NL. EQ. NDTOP) THEN 
P(J,JJ) = -UBCOND*FROW(L,KND,K) + P(J,JJ) 
ELSEIF(NL. EQ. NDTOP+1) THEN 
P( J, PERND- JJ+1 )=-UBCOND**FROW(L,KND,K)+P( J, PERND-JJ+1) 
ELSEIF((JJ. EQ. (BIND+1)). AND. (NL. EQ. NDTOT)) THEN 
C(J,1) = FROW(L,KND,K) + C(J,1) 
ELSEIF(NL. EQ. NDTOT) THEN 
P(J,JJ+1) = -UBCOND*FROW(L,KND,K) + P(J,JJ+1) 
6 UPPER ROW 
ELSEIF(NL. LT. NDTOP) THEN 
B(J,NL-1) = FROW(L,KND,K) + B(J,NL-1) 
c LOWER ROW 
ELSEIF(NL. LT. NDTOT) THEN 
C(J,NL-NDTOP-1) = FROW(L,KND,K) + C(J,NL-NDTOP-1) 
C ERROR 


Bick 
WRITE(*,*) NL,'ERROR - DUE TO INDEX GREATER THAN NODE TOTAL ' 
ENDIF 
C 
600 CONTINUE 
700 CONTINUE 
800 CONTINUE 
C 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(1,*) JJ 
DO 91, M = 1, NABC(JJ,2) 
WRITE(1,1020) (REAL(A(M,N)), N 
91 CONTINUE 
WRITE(1,*) ° ' 
DO 92, M = 1, NABC(JJ,2) 
WRITE(1,1020) (REAL(B(M,N)), N 
92 CONTINUE 
WRITE(1,*) ° ' 
DO 93, M = 1,: NABC(JJ,2) 
WRITE(1,1020) (REAL(C(M,N)), N 
93 CONTINUE 
WROTE Gis?) ae 
DO 94, M = 1, NABC(JJ,2) 
WRITE(1,1020) (REAL(P(M,N)), N 
94 CONT INUE 


1, NABC(JJ,1)) 


1, NABC(JJ,2)) 


1 eNGeGGII, 3) 


1, PERND) 


WRITE(1,*) ° ' 

WRITE(1,*) [ ' 

WRITE CIE. meee 
ENDIF 


IF(. NOT. ((CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 


CALL MARCH( JJ, IMX,NBMX ,NABC(JJ,1),NABC( JJ,2) ,NABC(JJ,3),1NORM, 
CSVEC) 
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CALL ZERO 
ENDIF 


CONTINUE 
LOAD A, B & P FOR I = BIND + 2 (BOTTOM) 


J=1 
ioe50, JD = 1, NCT(7) 
fe= NOC] ID) 
POO. kK =. 93 
ND(K) = LND(L,K) 
IF(CND(K).EQ.7) KND = K 
CONTINUE 


BOUNDARY CONDITION 
IFCNL. EQ. 6) THEN 
P(J,JJ+1) = -UBCOND*FROW(L,KND,K) + P(J,JJ+1) 
UPPER ROW 
ELSEIF((NL. GT. 1). AND. (NL. LT. 5)) THEN 
A(J,NL-1) = FROW(L,KND,K) + A(J,NL-1) 
LOWER ROW 
ELSEIF(NL. EQ. 7) THEN 
B(J,1) = FROWC(L,KND,K) + B(J,1) 
ERROR 
ELSE 
WRITE(*,*) NL,'ERROR - INDEX EQUAL TO 1 OR 5' 
ENDIF 


CONTINUE 
CONTINUE 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(1,*) TCALL+1 
DO 61, M = 1, NABC(TCALL+1,2) 
WRITE(1,1020) (REAL(A(M,N)), N 
CONTINUE 
WRITE(1,*) ' ' 
DO 62, M = 1, NABC(TCALL+1,2) 
WRITE(1,1020) (REAL(B(M,N)), N 
CONTINUE 
WRITE(1,*) ' ' 
DO 64, Mest, NABCCTCALLEL, 2) 
WRITE(1,1020) (REAL(P(M,N)), N 
CONTINUE 


iene pe TCALLEL 1) ) 


1, NABC( TCALL+1,2)) 


Is ASE 


RRUPEGIcs) 9) 5! 

WRIBE(1,*) ° ‘ 

WRITE(1,*) ‘ ' 
ENDIF 


WRITE(6,*) ' FINAL INVERSION' 

IF(. NOT. (( CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm'))) THEN 

CALL MARCH( JJ+1, IMX,NBMX , NABC( JJ+1,1) ,NABC( JJ+1,2) ,NABC( JJ+1,3), 
CINORM,SVEC) 
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ENDIE 


WRITE(2,*) . END END ' 


WRITE(6,*) 

WRITE(6,*) ‘MINIMUM AREA = ' ,MINAREA 

WRITE(6,*) ‘AT ROW ' sMINROW , ' AND ELEMENT NUMBER ' ,MINEL 
WRITE(6,*) 'MAXIMUM AREA = ° ,MAXAREA 

WRITE(6,*) 'AT ROW ' ,MAXROW, AND ELEMENT NUMBER ' ,MAXEL 
RATIO = MAXAREA/MINAREA 

WRITE(6,*) 'AREA RATIO = ',RATIO 


IF(RATIO. GT. 2.5) THEN 
WRITE(6,*) 'YOU SHOULD CONSIDER ABORTING THIS RUN AND ' 
WRITE(6,*) 'LOOKING AT THE MESH IN CURVE DIGITIZER ' 
WRITE(6,*) 'A BETTER METHOD MAY BE AVAILABLE ' 

ENDIF 

IF((CHAR3. EQ. 'D'). OR. (CHAR3. EQ. 'd')) THEN 
WRITE( 20,1080) 
WRITE( 20,1090) 
WRITE( 20, 1100) 

ENDIF 


FORMAT( 16,° OUL OR’ 3i6, Se7rus = 
FORMAT( I16,' OUT OF',I6,' CALLS ') 
FORMAT( 20( E14. 8,1X,E14. 8,1X)) 
FORMAT( 6X,'CALL COMPRS' ) 

FORMAT( 6X,'CALL NOBRDR’ ) 

FORMAT( 6X,'CALL PAGE(8.0,10.0)') 
FORMAT( 6X, 'CALL AREA2D(5.0,7.0)') 
FORMAT( 6X,'CALL FRAME' ) 

FORMAT( 6X,'CALL DONEPL' ) 

FORMAT( 6X,'STOP' ) 

FORMAT( 6X,'END') 


RETURN 
END 


SUBROUTINE FILL(I,LEL,FROW, ALPHA, BETA) 


COMPLEX F(3,3), FROW(100,3,3), ALPHA, BETA, A, B 

REAL AREA, MINAREA, MAXAREA 

REAL NDP( 200, 2) 

INTEGER I, J, K, L, LEL, LND(0:200,3), NDL(200,4), NCT(200) 
INTEGER MINROW, MINEL, MAXROW, MAXEL, M 

CHARACTER*1 CHAR2, CHAR3, CHAR4 

COMMON/BLK4/LND ,NDL,NCT 

COMMON /BLK5 /NDP 

COMMON /BLK7 /MINAREA , MINROW , MINEL, MAXAREA , MAXROW, MAXEL, AREA 
COMMON/BLK9/CHAR2, CHAR3, CHAR4 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(11,%*) ‘ROW NUMBER = ',I 

ENDIF 

DO 30 28) —s las 
IF((CHAR4. EQ. 'U'). OR. (CHAR4. EQ. 'u')) THEN 
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19 
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He 


ENDIF 
ENDIF 
CALL VARINT(J,F,A4,B,AREA, LND) 
ined. Gho2). AND, (J. Lf. LEL-1)) THEN 
IF( AREA. LT. MINAREA) THEN 
MINAREA = AREA 
MINROW = I 
MINE Ee — a 
ELSEIF( AREA. GT. MAXAREA) THEN 
MAXAREA = AREA 
MAXROW = I 
MAXEL = J 
ENDIF 
ENDIF 
DO 20, K = 1, 3 
WRITE(2,*) NDP(LND(J,K),1), NDP(LND(J,K),2) 
DO! WG he 
ROW vic Lye =e Gkeul) 
CONTINUE 
DEC CCHAR?2 EO moe OR CCHARZ EQ. i ))} THEN 
VWelin Gil 100 es kK RbAlLen(K,L)), & = 1,3) 
ENDIF 
CONTINUE 
DMG@cCHAR?. HO. Is), OR COHAR2. EO. a )) THEN 
WRITE(11,*) ° ' 
ENDIF 
WRITE(2,*) NDP(LND(J,1),1), NDP(LND(J,1),2) 
WRITE(2,*) ' 999990 999990 ' 
DISSPLA PROGRAM GENERATION 
IF((CHAR3. EQ. 'D'). OR. (CHAR3. EQ. 'd’)) THEN 
WRITE( 20,1010) NDP(LND(J,1),1), NDP(LND(J,1),2) 
WRITE( 20,1020) NDP(LND(J,2),1), NDP(LND(J,2),2) 
WRITE( 20,1020) NDP(LND(J,3),1), NDP(LND(J,3),2) 
WRITE( 20,1020) NDPC(CLND(J,1),1), NDPCLND(J,1),2) 
ENDIF 


CONTINUE 


HOnRMANG@ ie 15,2X),3X,3(CF8.5,2X)) 


BORMAR@OM. CALL SERTPT(' ,F8.5,',',F8.5,°)') 
FORMAT( 6X, ‘CALL CONNPT(',F8.5,',',F8.5,')') 
RETURN 

END 
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SUBROUTINE ZERO 

ZERO A, B, C, AND P MATRICES 

COMPLEX A(50,50), B(50,50), C(50,50), P(50,100) 
INTEGER J, K, L 


COMMON/BLK8/A,B,C,P 


DO 50, J = 1, 50 
DO 400k = sO 


A(J,K)} = CMPLX(0. 0,0. 0) 
B(J,K) = CMPLX(0.0,0. 0) 
C(J,K) = CMPLX(0.0,0.0) 
CONTINUE 
DO 60, L= 1, 100 
P(J,L) = CMPLX(0.0,0. 0) 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE BNDC(I,BCOND,E0O,SURBC,CHAR1,ALPHA, BETA, LINE,MODE, 


CXORIGIN, YORIGIN, CHAR4) 


BOUNDARY CONDITION FILL FOR A PLANE WAVE 
EO = FIELD STRENGTH 
KO = WAVE NUMBER (2*PI/WAVELENGTH) 


BOUNDARY CONDITION FILL FOR A CYLINDRICAL BOUNDARY CONDITION 
EO = FIELD STRENGTH 
MODE = MODE NUMBER FOR CYLINDRICAL BOUNDARY CONDITIONS 


COMPLEX SURBC(100), VALUE(50), BCOND(100), ALPHA, BETA, LINE(50) 
COMPLEX K, RAK 

REAL NDP(200,2), EO, KR, KI, XORIGIN, YORIGIN, RA, RB 

INTEGER I, J, NDTOP, NDBOT, NDTOT, NODES(200), PERND, BIND, MODE 
CHARACTER*1 CHAR1, CHAR4 


COMMON/BLK2/PERND, BIND 
COMMON/BLK3/NODES 
COMMON/BLK5/NDP 


IF((CHAR1. EQ. 'P'). OR. (CHAR1. EQ. 'p')) THEN 
IF((CHAR4. EQ. 'U'). OR. (CHAR4. EQ. 'u')) THEN 


KR = REALCCSQRT( BETA) ) 

KI = ABS(CAIMAG( CSQRT( BETA) ) ) 
ELSE 

KR = 1.0 

KI = 0.0 
ENDIF 


CHECK ON WHETHER WE ARE USING ER = ALPHA OR BETA 


IFC 1. nO] 1) THEN 
BCOND( 1) = EO*“EXPCKI*NDP(1,2))*CMPLX(COSCKR*NDP(1,2)), 
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e SIN(KR*NDP(1,2))) 
BCOND(PERND) = EO*EXP(KI*NDP(3, 2) )*CMPLX(COS(KR*NDP(3,2)), 
C SIN( KR*NDP(3,2))) 
BCOND(2) = EO*EXP(KI*NDP(7,2))*CMPLX(COS(KR*NDP(7,2)), 
C SIN(KR*NDP(7,2))) 
VALUE(1) = EO**EXP( KI*NDP( 2,2) )*CMPLX(COS(KR*NDP(2,2)), 
a SIN(KR*NDP(2,2))) 
WRITE(10,*) BCOND(1) 
WRITE(10,1000) VALUE(1) 
WRITE(10,*) ' ' 
SURECOME— VALUECL) 
ELSEIF(1I. LE. BIND) THEN 
NDTOP = NODES(I) + NODES(PERND-I+2) - 1 
NDBOT = NODES(I+1) + NODES(PERND-I+1) - 1 
NDTOT = NDTOP + NDBOT 
BCOND( PERND-I+1) = EO*EXP(KI*NDP(NDTOP+1, 2) )*CMPLX(COS(KR* 
ce NDP(NDTOP+1 ,2)),SIN( KR*NDP( NDTOP+1, 2) )) 
BCOND(I+1) = EO*EXP(KI*NDP(NDTOT, 2) )**CMPLX(COS(KR*NDP(NDTOT, 
C 2)),SIN(KR*NDP( NDTOT, 2))) 
WRITE(10,*) BCOND( PERND-I+2) 
DO 10, J = 2, NDTOP-1 
VALUE(J) = EO*EXP(KI*NDP(J,2))*CMPLX(COS(KR*NDP(J,2)), 
C SIN(KR*NDP(J,2))) 
IF(J.EQ. 2) THEN 
SURBC(PERND-I+2) = VALUE(2) 
ELSEIF(J. EQ. (NDTOP-1)) THEN 
SURBC(I) = VALUE(NDTOP-1) 
ENDIF 
CONTINUE 
LINE(I) = VALUE( (NDTOP+1)/2) 
WRITE(10,1000) (VALUE(J), J = 2, NDTOP-1) 
WRITE(10,**) BCOND(1I) 
WRITE(10,*) | ' 


ELSEIF(1I. EQ. BIND+1) THEN 

BCOND(I+1) = EO*EXP(KI*NDP( 6,2) )*CMPLX(COS(KR*NDP(6,2)), 
@ SIN(KR*NDP(6,2))) 

VALUE(2) = EO*EXP(KI*NDP(2,2))*CMPLX(COS(KR*NDP(2,2)), 
C SIN(KR**NDP(2,2))) 

VALUE(3) = EO*EXP(KI*NDP(3,2))*CMPLX(COS(KR*NDP(3,2)), 
CG SINC KR*NDP(3,2))) 

VALUE(4) = EO**EXP(KI*NDP( 4,2) )*CMPLX(COS(KR*NDP(4,2)), 
G SIN(KR*NDP(4,2))) 

WRITE(10,*) BCOND(I+2) 

WRITE( 10,1000) VALUE(2), VALUE(3), VALUE(4) 

SURBC(BIND+3) = VALUE(2) 

SURBC(BIND+1) = VALUE(4) 

LINE(I) = VALUE(3) 

WRITE(10,*) BCOND(I) 

WRITE(10,*) ° ' 

VALUE(2) = EO*EXP(KI*NDP(7,2))*CMPLX(COS(KR*NDP(7,2)), 
é SIN(KR*NDP(7,2))) 

WRITE(10,1000) VALUE(2) 

WRITE(10,**) BCOND(I+1) 

SURBC(BIND+2) = VALUE(2) 
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ENDIF 
ELSEIF((CHAR1. EQ. 'C'). OR. (CHAR1. EQ. 'c')) THEN 


IF(I. EQ. 1) THEN 
RA = SQRT((NDP(2,1)-XORIGIN)***2+(NDP( 2,2)-YORIGIN)**2) 
RB = SORT((NDP(1, 1)-XORIGIN)**2+(NDP( 1,2)-YORIGIN)**2) 
WRITE(**,*) BETA 
K = CSQRT( BETA) 
RAK = RA*K 
CALL CYLBC(RA,RB,RAK,NDP(1,1),NDP(1,2),XORIGIN,YORIGIN,E0, 
BCOND( 1) ,SURBC(1),K) 
CALL CYLBC(RA,RB,RAK,NDP(3,1),NDP(3,2),XORIGIN, YORIGIN,EO, 
BCOND( PERND) , SURBC( PERND) , K) 
CALL CYLBC(RA,RB,RAK,NDP(7,1),NDP(7,2),XORIGIN, YORIGIN,EO, 
BCOND( 2) , SURBC(2) ,K) 
WRITE(10,%*) BCOND(1), SURBC(1) 
WRITE(10,*) ' 

ELSEIF(I. LE. BIND) THEN 


NDTOP = NODES(I) + NODES(PERND-I+2) - 1 
NDBOT = NODESCI+1) + NODES CPERND ier al 
NDTOT = NDTOP + NDBOT 


CALL CYLBC(RA,RB,RAK,NDP(NDTOP+1,1),NDP(NDTOP+1,2),XORIGIN, 
YORIGIN , EO , BCOND( PERND-I+1) , SURBC( PERND-I+1) ,K) 
CALL CYLBC(RA,RB,RAK,NDP(NDTOT, 1),NDP(NDTOT,2),XORIGIN, 
YORIGIN, EO, BCOND( I+1) , SURBC( I+1) ,K) 
WRITE(10,%*) BCOND( PERND-I+2), SURBC( PERND-I+2) 
WRITE(10,**) BCOND(I), SURBC(1) 
WRITE(10,*) ' ' 
ELSEIF(I.EQ. BIND+1) THEN 
WRITE(10,%*) BCOND(I+2), SURBC(I+2) 
WRITE(10,*) BCOND(I), SURBC(I) 
WRITE(10,%*) ' 
CALL CYLBC(RA,RB,RAK,NDP(6,1),NDP(6,2),XORIGIN, 
YORIGIN, EO, BCOND( I+1) , SURBC( I+1),K) 
WRITE( 10 ,**) BCOND(I+1), SURBC(I+1) 
WRITE(10,**) 
ENDIF 


ELSE 
RETURN 
ENDIF 
FORMAT( 1X,50(E14. 8,2X,E14. 8,2X)) 
RETURN 
END 
SUBROUTINE CYLBC(RA,RB,RAK,X,Y,XORIGIN, YORIGIN, EO, BC,PSI,K) 
COMPLEX SQRTM1, JORB, J1RB, JORAK, J1RAK, HORA, H1RA, HORB, H1RB 
COMPLEX DELTAN, AN, PSI, JORA, J1RA, K, RAK, BC 
REAL PI, XORIGIN, YORIGIN, X, Y, RA, RB, EO, PHI 


PT = 4. O*ATAN(C 1.0) 
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SQRTM1 = CMPLX(0. 0,1. 0) 

PHI = ATAN2(X - XORIGIN, Y ~ YORIGIN) 
CALL BES1(CMPLX(RA,0. 0) , JORA, J1RA) 
CALL BESI(CMPLXCRB,O.0) ,JORB,J1RB) 
CALL BES1(RAK, JORAK, J1RAK) 

CALL HAN1(CMPLX(RA,O. 0) ,HORA,H1RA) 
CALL HAN1(CMPLX(RB,0. 0) ,HORB,H1RB) 


DELTAN = J1RB*(J1RAK*( HORA-H1RA/RA) -K*( JORAK-J1RAK/RAK)*H1RA) - 


CH1RB*( J1RAK*( JORA-J1RA/RA) -K*( JORAK-J1RAK/RAK)*J1RA) 
AN = -2. 0*SQRTM1/(PI*RA*DELTAN) 

BC = EO*COS( PHI) 

PSI = AN*J1RAK*BC 


RETURN 
END 


‘oe. 


SUBROUTINE BESI(Z,J0,J1) 
Computing Bessel Functions for n = 0, 1 with 
CABS(Z) .LE. 6 and Hankel'’s Asymptotic Formula for 


GABSCZ) .GT. 6. 
Written 11/5/87 by M.A. Morgan 


Gi) C? CrGicc] 


INTEGER M,M2 

REAL C(34) ,DM,F(34) ,GO,P(34) ,Pi,P2 

COMPLEX Z,Z2,23,24,J0,J1,AM,CL,PO,P1,Q0,Q1,C0,C1,S80,§ 
Pi=3. 1415927 

P2=2.0/PI 

IF(CABS(Z). LE. 6.0) THEN 


Utilizing the Direct Power Series Method 


i lp fa | 


GO= 1.781072 
Z2=0.5*Z 
CL=CLOG(G0*Z2) 


Comoltimec hm) — m! and P(m) = 1+ 1/2 + 1/3 + .. 


eo Nis <p Ya a 


F(1)=1.0 
P(1)=1.0 
DO 11 M=2,34 
F(M)=M*F(M-1) 
P(M)=P(M-1)+1. 0/M 
1 CONTINUE 


Computing Power Series Coefficients 


qaqa 


DM=-1.0 
DO 22 M=1,34 
CCM)=DM/C FCM) *F(M) ) 
DM=-DM 
22 CONTINUE 
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Complex Argument Z. Direct Power Series Method for 


1 


aly 1 


5s 


C) C76) C2 Cao) 


Co C37@) 


Computing JO and Jl 


JO=(1. ,0. ) 
J1=(0. ,0. ) 
M=0 
M=M+1 
M2=2°M 
AM=C(M)**( Z2***M2) 
JO=JO+AM 
J1=J1-M*AM 
IF((CABS( AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 33 
J1=31/22 
return 

ELSE 


Hankel’ Asymptotic Formula (Abram. & Stegun p. 364) 


Z2=2*7Z 
73=7570 
Z4=2723 
PO=1. 0-. 0703125/Z22+. 1121521/24 
QO=-. 125/2+. 0732422/23 
P1=1. 0+. 1171875/Z2-. 1441956/24 
Q1=. 375/Z-. 10253906/Z3 
CO=CCOS(Z-. 25**PT) 
SO=CSIN(Z-. 25*PI) 
C1=CCOS(Z-. 75**PT) 
S1=CSIN(Z-. 75*PI) 
AM=CSQRT(P2/Z) 
JO=AM**( PO*C0-Q0*S0) 
J1=AM**(P1*C1-Q1*S1) 

ENDIF 

RETURN 

END 


SUBROUTINE HAN1(Z,HO,H1) 


Computing Hankel Functions for n = 0, 1 with 

Complex Argument, Z. Direct Power Series Method for 
CABS(Z) .LE. 5 and Hankel's Asymptotic Formula for 
CABS CA). Glo. Written 11/6/87 by M.A. Morgan 


INTEGER M,M2 

REAL C(34),DM,F(34) ,GO,P(34),Pi,P2 

COMPLEX Z,22,23,24,J0,J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 
COMPLEX E0,E1,X0,X1,HO0,H1, j 

PI=3. 1415927 

P2=2. 0/Pi 

j=(0. ,1.) 

IF(CABS(Z). LE. 5.0) THEN 


Direct Power Series Method 


GO= 1. 78072 
Z2=0. 5*Z 
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33 


2 © 


44 


Computing F(m) =m! 


Computing Power Series Coefficients 


CL=CLOG(G0*Z2) 


al) — io 

ea j=1, © 

DO 11 M=2,34 
F(M)=M*F(M-1) 


ands ecamje= let 1/2 +-1/3 +... 
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CONTINUE 


DM=-1. 0 
DO 22 M=1,34 


se SS) 


DM=-DM 
CONTINUE 


Computing JO and Jl 


JO=(1. ,0.) 
J1=(0. ,0. ) 

M=0 

M=N+1 

M2=2"M 

AM=C(M)**( Z2***M2) 
JO=JO+AM 
J1=J1-M*AM 


sg ar yan 


PeCeABoCAM Gil i. 0OE-10).AND)(M.bT.34)) GO TO 33 


J1=J51/2Z2 


Computing YO and Yl 


ELSE 


Hankel' 


M=0 

YO=CL* JO 
Y1=Z2*CL*J1-0. 5*J0 
M=M+1 

M2=2*M 
AM=C(M)*P(M)**( Z2****M2 ) 
YO=YO-AM 

Y1=Y1+M*AM 


IF((CABS( AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 44 


YO=P2**Y0 
Y1=P2*Y1/Z2 
HO=JO-j*YO 
H1=J1-j*Y1 
RETURN 


Asymptotic Formula (Abram. & Stegun p. 364 


DDN 
73-2 
Z4=Z2*Z3 


PO=1. O-. 0703125/Z2+. 1121521/24 


itl 


ClrC 2G) GC) 


C2C2-€2 €2 C200 


UNE 


i 


QO=-. 125/Z+. 0732422/23 

P1=1. O+. 1171875/Z2-. 1441956/Z4 
Q1=. 375/Z-. 10253906/Z3 
XO=(Z-. 25*PT) 
X1=(Z-. 75*PI) 
EO=CEXP( - j**X0) 
E1=CEXP( - j**X1) 
AM=CSQRT( P2/Z) 
HO=AM**( PO- 4**Q0)*E0 
H1=AM*(P1-4*Q1)*E1 

ENDIF 

RETURN 

END 


SUBROUTINE MARCH( I, IMX,NBMX,NA,NB,NC,INORM,SVEC) 
THIS ROUTINE PERFORMS THE RICCATI TRANSFORM FIRST 


SWEEP, GENERATING AND STORING ON DISK 1 RMAT 
AND SVEC (FOR EACH MODE) AT EACH FORWARD STEP. 


COMPLEX RMAT(50,50),SVEC(50,100) 

COMPLEX A(50,50),B(50,50),C(50,50),P(50, 100) 
COMPLEX D(50),SUM,DET 

REAL COND, DMAG 

INTEGER I,J,K,L,NA,NB,NC,PERND,BIND, IMX, INORM 
INTEGER NBMX 

CHARACTER*1 CHAR2, CHAR3, CHAR4 


COMMON /BLK2/PERND , BIND 
COMMON /BLK8/A,B,C,P 
COMMON/BLK9/CHAR2, CHAR3, CHAR4 


LOADING THEN INVERTING (B+A*RMAT) 
USING MINIMUM MEMORY SINGLE MATRIX TECHNIQUE 
DATE 1/29/80 FOR THIS CHANGE 


SKIPPING FIRST A*R (WHEN A = 0, ZERO R MATRIX) 
IF(I. EQ. 1) THEN 

DO 10, J = 1, NB 

WO, KS fh, ING 
WP RMATC I,K) = (0. ,0.) 

CONTINUE 
RMAT = A*R 
ELSE 


IF(( CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(19,*) ' OLD R MATRIX' 
er ia, Je i. Ss 


WRITE( 19,1000) CREALCRMAT(J,K)), K = 1, 5) 


CONTINUE 

WRITE(19,*) ' ' 
WRITE(19,*) ' A MATRIX’ 
DO 1223) =i 


WRITE( 19,1000) CREALCA(J,K)), K = 1, 5) 


CONTINUE 
WRITE(19,*) ' ' 
ENDIF 


lel 


20 


ro 


14 


iS 


ap las <p Ta > 1 > | 


DOmSOeek = 1, NB 


CONTINUE 
HO 30, J ="1, NB 
RMATC I,K) = DCJ) 
CONTINUE 


IF( (CHAR2. EQ. 'I'). OR. (CHAR2.EQ. 'i')) THEN 
WRITE(19,*) ' NEW R MATRIX’ 
Doeow = 1, 5 


WhhPe¢19. J0G0) (REAECKRMAT(CJ,K)), K = 1, 5) 


CONTINUE 
WRITE(19,*) ° ' 
ENDIF 


END IE 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(19,*) ' B MATRIX' 
Omics = ls 5 
WRITE( 19,1000) CREALCB(J,K)), K = 1, 5) 
CONTINUE 
WRITE(19,*) ' ' 
ENDIF 


B + RMAT 
DO 40, J = 1, NB 
DO 40, K = 1, NB 
RMAT(J,K) = RMAT(J,K) + B(J,K) 
CONTINUE 


IF((CHAR2. EQ. 'I'). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(19,*) ' NEWEST R MATRIX' 
Domo. 7 = ene 


WRITE(9,1000) (REALC(CRMAT(J,K)), K = 1, NB) 
WRITE( 19,1000) (REALCRMAT(J,K)), K = 1, NB) 


CONTINUE 
WRITE(9,*) ° ' 
WRITE(19,*) ' ' 
ENDIF 
INVERTING THE MATRIX (B + A*R) 
CALL CSMINV(RMAT ,NBMX,NB ,DET,COND , INORM) 


DMAG = CABS(DET) 
WRITE(*,*) I,NBMX,NB,DMAG,COND 


COMPUTING THE NEW S-VECTORS 


SKIPPING FIRST A*S (A = 0) 
PC Ee) thin 
CONTINUE 
ELSE 
DOR Oj — I PERND 


eS 


DO 60, J = 1, NB 
D(J) = (0.0,0.0) 
DO 60, L = 1, NA 
D(J) = D(J) + ACJ,L)*SVEC(L,K) 
60 CONTINUE 
DO 70, J = 1, NB 
P(J,K) = PCs, 
70 CONTINUE 
ENDIF 
C FINAL SVEC MULTIPLICATION 
DO 100, K = 1, PERND 
DO 90, 0 = 1, Ne 
D(J) = (0.0,0.0) 
DONOGmn iy eel 
D(J) = D(J) + RMAT(J,L)*P(L,K) 
90 CONTINUE 
DO 100, J = 1, NB 
SVEC(J,K) = D(J) 


100 CONTINUE 
C 
C SIGRING T+1 SVEGSGN Disha? 


DO 110° J = 1 NE 
WRITE(7) (SVECCJ,K)5 °K = 1) PERND) 
10 CONTINUE 


C 
¢ FINAL RMAT MULTIPLICATION 
IF(I.EQ. IMX) RETURN 
DO 130, J = 1, NB 
DO 120, K = 1, NC 
D(K) = (0.0,0.0) 
Nh i, NE 
D(K) = D(K) - RMAT(J,L)*C(L,K) 
120 CONTINUE 


DO 130, K = 1, NC 
RMAT(J,K) = D(K) 
130 CONTINUE 


a 


2 STORING I+1-RMAT ON DISK 7 

DO 140, J = 1, NB 

WRITE(7) (RMAT(J,K), K = 1, NC) 
140 CONTINUE 


C 
1000 FORMAT(10(E9. 3,1X)) 
C 
RETURN 
END 
E 
E 
SUBROUTINE SWEEP( IMX,NABC,SURBC, LINE,CHAR1,U,BCOND,PSI,ANS, CHARS) 
€ THIS ROUTINE PERFORMS THE RICCATI TRANSFORM BACKSWEEP 
c FROM I=IMX TO I=1, RECALLING RMAT AND 
c SVEC FROM DISK 7 AT EACH BACKSTEP TO FORM 
C THE NODE VECTORS PSI, THEN STORING THESE ON 
c DISK 8 FOR EACH APPLIED DIRICHLET B.C. 


COMPLEX RMAT(50,100),PSI(50,100),SVEC(50,100) ,D(100),TEMP 
COMPLEX SURBC(100),ANS( 100) ,LINE(50),ANSB(50),U( 100,100) 
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COMPLEX BCOND( 100) 

REAL ERRORP,ERRORD,TERRN, TERRD,ATSERR 
iViEGRReI J j,k NABCC100,3) ,1MX,PERND,BIND 
CHARACTER*1 CHAR1, CHARS5 


COMMON/BLK2/PERND, BIND 


IF( (CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm')) THEN 
RETURN 
ENDIF 


INITIAL DISK READ AT IMX (R = 0, NOT WRITTEN, => S IS READ FIRST) 
WRITE(*,*) ' ' 
hoeo0. I = IMX, 1, -1 
WRITE(*,1050) (IMX-I+1), IMX 
IF(I.EQ. IMX) THEN 
Domlee Je— NABC(T.2)5 1, =1 
BACKSPACE 7 
READ(7) G@PSI(J,K), K = 1, PERND) 
BACKSPACE 7 
CONTINUE 
DOw20 59) — 1 NABCCI.2) 
DO 15, K = 1, PERND 
Vente ky = PSI(J.K) 
CONTINUE 
CONTINUE 
SUBSEQUENT DISK READS 
ELSE 
READ R MATRIX 
DO 30, J = NABC(I,2), 1, -1 
BACKSPACE 7 
READ(7) (RMAT(J,K), K = 1, NABC(I,3)) 
BACKSPACE 7 
CONTINUE 
READ S VECTOR 
PO 40,00 = NABC( 1,2), 1, -1 
BACKSPACE 7 
READ(7) (SVEC(J,K), K = 1, PERND) 
BACKSPACE 7 
CONTINUE 
MULTIPLY RMAT = RMAT*PSI 
DO 60, J = 1, NABC(I1,2) 
DO 50, K = 1, PERND 
CK) = Col, (ey ole) 
DOwsOr Ly = IeeNAbecl 3) 
D(K) = DCK) + RMAT(J,L)*PSICL,K) 
CONTINUE 
BOmoO). K —eles PERND 
RuAamCIeR) = DOK) 
CONTINUE 
PSI = RMAT + SVEC 
DOO me— ele NABCCI 2) 
DO 70, K = 1, PERND 
PSI(J,K) = RMAT(J,K) + SVEC(J,K) 
CONTINUE 
ANSB(1) = (0.0,0.0) 


Ds 


80 


90 


94 
98 


100 


110 


DO 80, J = 1 NABCGI 2) 
DO 75,08 eee ND 

IFC (J. EQUNABCCL 2 Ok lee) aun 
UCT SK) = PsiCsa 

ELSEIFC J. EQ: 1) THEN 
UC 2" 1MA-1,K) = ees jiae 

ELSETFCJ. EQ. (NABC( 1,2) 41) 2) etn 
ANSBCI) = ANSBCI) + PSI(J,K)*BCOND(K) 


ENDIF 
CONTINUE 
CONTINUE 
ENDIF 
CONTINUE 
WRITE(8,*) ° ' 


DO 98, J = 1, BERND 
ANS(J) = (0.0,0. 0) 
DO 94, L = 1, PERND 
ANS(J) = ANS(J) + UCJ,L)*BCOND(L) 


CONTINUE 

CONTINUE 

TERRN = 0.0 

TERRD = 0.0 

DO 100, I = 1, PERND 
ERRORP = (CABS(ANS(I) - SURBC(I)))**2 
ERRORD = (CABS(SURBC(I)))7**2 


TERRN = TERRN + ERRORP 
TERRD = TERRD + ERRORD 
WRITE( 13,1020) I, BCOND(I), SURBC(I), ANS(I), ERRORP 
CONTINUE 
WRITE(13,*) | ' 
ATSERR = (TERRN/TERRD)***0. 5 
WRITE( 13,1030) ATSERR 
WRITE(*,*) "RMS ERROR (FOR THE PERIMETER) = ' ,ATSERR 


IF((CHAR1. EQ. "C'). OR. (CHARI. EQ. 'c')) THEN 
RETURN 
ELSE 
TERRN 
TERRD 
DOr 10% 
ERRORP = (CABS(ANSB(I) - LINE(I)))**2 
ERRORD = (CABS(LINE(1)))**2 
TERRN = TERRN + ERRORP 
TERRD = TERRD + ERRORD 
WRITE( 13,1025) (I-1), LINE(I), ANSB(I), ERRORP 
CONTINUE 
WRELECI 30a 
ATSERR = (TERRN/TERRD)**0. 5 
WRITE( 13,1040) ATSERR 
WRITE(*,*) 'RMS ERRROR (FOR BISECTION SEGMENT) = ',ATSERR 
ENDIF 


0.0 
0.0 
le= 2, BENDA 


DO 120, I = 1, PERND 
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—-o-aremnoung — ae 


CI1GarG) 
Oo 


Cy Ciena C7 C2 C2°C) C2 C2 C9") 


CT W Ph 


WeinesOr O60) (UCT .J), J.= 1, PERND) 
CONTINUE 


mORMATC IN, 13,2X,3(614. 8,1X,B14. 8,3X),F10. 6) 
FORMAT(1X,13,2X,4(E14. 8,2X) ,F10. 6) 

FORMAT(1X,' RMS ERROR (FOR THE PERIMETER) = 12a) 
FORMAT(1X,' RMS ERROR (FOR BISECTION SEGMENT) = ' /,F12. 6) 
FORMAT(1X,13,' TOTAL BACKSWEEP ROWS OUT OF ',13,' COMPLETED’ ) 
FORMAT( 1X, 100(E8. 2,E8. 2,2X)) 


RETURN 
END 


SUBROUTINE CSMINV(CA,NDIM,N,DETERM, COND, INORM) 


INORM - FLAG TO NORMALIZE COLUMNS AND ROWS OF MATRIX A 
MATRIX NORMALIZATION BY M.A. MORGAN 
APRIL 24,1978 


A e Math x10 INVERT =I NEU OUTPUT 
NDIM = SNE UT 
N 9 SNE OT 
WRTERM - DETERMINATE OF A SOU IE Ut: 
COND - CONDITION NUMBER OF A secUTPUT 


INORM - INTEGER NORMALIZATION FLAG - INPUT 


COMPLEX A(50,50),PIVOT(50) ,AMAX,T, SWAP, DETERM,U 

INTEGER I,J,K,L,IPIVOT(50) , INDEX(50,2) , IROW, ICOLUM,L1, JROW 
INTEGER JCOLUM,N, INORM 

REAL TEMP, ALPHA(50) ,COL(50) , ROW(50) , AJK, SUMAXA, SUMROW , SUMAXI 


IF(NDIM. GT. 50) THEN 
WRITE(**,*) ' ERROR IN INVERTION CALL... DIMENSION > 50 ' 
STOP 

ENDIF 

IF(N. GT. NDIM) THEN 
WRITE(*,*) ' ERROR IN INVERTION CALL... N > MAX DIM. ' 


STOP 
ENDIF 
IF(INORM.NE.1) GO TO 7 
DO 3 K=1,N 
eonCK = 10e0 
POmee— 1.N 


AJK = CABS(A(J,K)) 
TRCAIK.GT. COL(K)) COLCK) = AJK 
CONTINUE 
DO2J=1,N 
Ge — AGIaK) (COLE) 
CONTINUE 
ROW NORMALIZING 
Dine = 1 
ROW(J) 
DO 4K 


WoW A 


Loe 


by, 


AJK = CABS(ACJ,K)) 
IFC AJK. GT. ROWCJ)) ROW(J) = AJK 
+ CONTINUE 
DO5 K= 


en 
ACJ,K) = 


A( J,K)/ROW(J) 
CONTINUE 
CONTINUE 
DETERM = CMPLX(1.0,0.0) 
SUMAXA = 0.0 
DO 20 J=1,N 
ALPHA(J) = 
SUMROW = 0. 
DO 10 I=1, N 
ALPHA( J) = ALPHA(J) + A(J,I)* CONJG(A(J,I)) 
10 SUMROW = SUMROW + CABS(A(J,I)) 
ALPHA(J) = SQRT(ALPHA(J) ) 
IF(SUMROW. GT. SUMAXA) SUMAXA = SUMROW 
20 IPIVOT(J) = 0 
DO 600 I =1, N 


SON 


G70 
0 


C 
C 
AMAX = CMPLX(0. 0,0. 0) 
DO 105 J=1, N 
IF (IPIVOT(J)-1) 60, 105, 60 
60 DO 100 K = 1, N 
IF (IPIVOT(K)-1) 80, 100, 740 
80 TEMP = AMAX**CONJG(AMAX) - A(J,K)*CONJG(A(J,K)) 
IF(TEMP)85,85,100 
85 IROW = J 
ICOLUM = K 
AMAX = A(J,K) 
100 CONTINUE 
105 CONTINUE 
IPIVOT(ICOLUM) = IPIVOTCICOLUM) + 1 
G 
E 


IF (IROW-ICOLUM) 140, 260, 140 
140 DETERM = -DETERM 
DO 200 L=1, N 
SWAP = A(IROW,L) 
AC IROW,L) = AC ICOLUM,L) 
200 ACICOLUM,L) = SWAP 
SWAP = ALPHA(IROW) 
ALPHA(IROW) = ALPHA(ICOLUM) 
ALPHA(ICOLUM) = SWAP 
260 INDEX(I,1) = IROW 
INDEX(I,2) = ICOLUM 
PIVOT(I) = AC ICOLUM, ICOLUM) 
U = PIVOT(1) 
DETERM = DETERM*U 
DETERM = DETERM/ALPHA( ICOLUM) 
TEMP = PIVOT(1)**CONJG( PIVOT(I)) 
IF( TEMP) 330,720,330 


530 AC ICOLUM,ICOLUM) = CMPLX(1. 0,0. 0) 
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Diese OsL = 1,N 


U = PIVOT(I) 

350 AC ICOLUM,L) = ACICOLUM,L)/U 

6 

c 

380 DO 550 Li = 1, N 
IF(L1-ICOLUM) 400, 550, 400 

400 T = ACL1, ICOLUM) 
A(L1,ICOLUM) = CMPLX(0.0,0.0) 
DO 450 L= 1, N 

U = AC ICOLUM,L) 

450 AGiie l= ACI.) = UFT 

550 CONTINUE 

600 CONTINUE 

c 

e 


foe DO 710 I = 1, N 
ih SN So a 
IF (INDEX(L,1) - INDEX(L,2)) 630, 710, 630 
630 JROW = INDEX(L,1) 
JCOLUM = INDEX(L,2) 
Demyooek = 1) N 
SWAP = A(K,JROW) 
A(K,JROW) = AC(K,JCOLUM) 
A(K,JCOLUM) = SWAP 
705 CONTINUE 
710 CONTINUE 
SUMAXI = 0.0 


DO 910 I=1, N 
SUMROW = 0.0 
DO 900 J=1, N 
900 SUMROW = SUMROW + CABS(A(I,J)) 


IF(SUMROW. GT. SUMAXI) SUMAXI = SUMROW 
910 CONTINUE 
COND = SUMAXA*SUMAXI 
IF(INORM. NE.1) GO TO 955 
DO 950 K=1, N 
DO 950 J=1, N 


950 AC(J,K) = ACJ,K)/CROW(K)*COL(J)) 
955 CONTINUE 
RETURN 


720 WRITE(*,730) 
730 FORMAT( ' MATRIX IS SINGULAR’) 
740 RETURN 

END 


SUEROUTINE SAVEC BCOND, ANS ,U,OFFSET,PERND,CHARS ,DPER,K,XORG, YORG, 
CNRES , MRES , LBIAS ,GBIAS , MAXD) 


THIS SUBROUTINE SAVES THE ESSENCE OF THE FINITE ELEMENT PROBLEM 
NeeA DATA PILE CALLED “ F3.DAT ". tHIS DATA IS NECESSARY TO 
SOLVE THE FIELD FEEDBACK FORMULATION, (F3). 


eo il Ge ba SP a i J 


COMPLEX BCOND(100),ANS(100),U(100,100) 
REAL OFFSET,MESH(0: 200,5),DPER,K,XORG, YORG,MRES , MAXD 
INTEGER PERND,I,J,NRES,LBIAS,GBIAS 


Ds 


CHARACTER*1 CHAR5 
COMMON/BLK1/MESH 


IF( (CHARS. EQ. 'M'). OR. (CHARS. EQ. 'm')) THEN 
RETURN 
ENDIF 


WRITE(40,*) PERND 
WRITE(40,*) OFFSET 
WRITE(40,**) DPER 
WRITE(40,*) K 
WRITE(40,%*) XORG 


WRITE(40,%) YORG 
WRITE(40,%) NRES 
WRITE(40,**) MRES 
WRITE(40,*) LBIAS 
WRITE(40,*) GBIAS 
WRITE(40,*) MAXD 


DO! 10ye t= ee 
DO 10, J = 1, PERND 
WRITE(40,%*) MESH(J,1) 
CONTINUE 


DO. 20, 9 — 1, Eee 
WRITE(40,*) BCOND(J) 
CONTINUE 


DO 30, J = 1, PERND 
WRITE(40,*) ANS(J) 
CONTINUE 


DO 40, I = 1, PERND 
DO 40, J = 1, PERND 
WRITE(40,*) U(I,J) 
CONTINUE 


RETURN 
END 
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mao 


Cc) Cr crc? G22 


APPENDIX D. VARINT CONVERGENCE PROGRAM 


TEST OF VARINT CONVERGENCE 


COMPLEX F(3,3), ALPHA, BETA, EXACT, CENTER, LEFT, RIGHT, TOP 
COMPLEX BOTTOM, SUM, CALC 

REAL X(3), Y(3), D, AREA, KR, PI, ERROR 

INTEGER I 

OPEN (UNIT = 1, FILE = 'C: MSFORT TEST. DAT', STATUS = ‘UNKNOWN’ ) 
BeHA = (1.0,0.0) 

BemA = (1.0,0.0) 

PI = 4. O*ATAN(1. 0) 


pO 5, I=1, 100 
D = 2. O*PI*“FLOAT(I)/100. 0 
KR = REAL(CSQRT( BETA) ) 


ft 


Cy eS Oe 
Gh = O40 
2) Se 
Me? = 040 
C2) = 010 
Y¥(3) =D 


CALL VARINT(X,Y,F,ALPHA, BETA, AREA) 


EXACT = CMPLX(COS(KR*0. 0) ,SINCKR*0. 0) ) 
CENTER = 4. 0*F(1,1) 
RIGHT = -2. 0*F(1,2)*CMPLX( COSC KR*Y(2)) ,SINCKR*Y(2))) 
Meee Ori) CE NC COSC KRY (3 0, oIN(C KR*Y(3))) 
ieee 2 Oe hee) OME ma eos Kkw1C2) ) jolIN(KR*Y(2))) 
BOGTOM = —=2- 0 hale 2) “CHP ENCeOS( -KR*Y(3)) ,SIN( -KR*Y(3))) 
sUhe= LOPS BOTTOM 2 LEFT + RIGHT 
CALC = SUM/CENTER 
ERROR = (CALC-EXACT) /EXACT 
WRITE(1,1000) I,D,EXACT,CALC,ERROR 

CONTINUE 


CLOSE( 1) 

FORMAT(1X,13,1X,F8. 5,1X,2(F8. 5, 1X,F8.5,1X) ,E13. 6) 

STOP 

END 

SUBROUTINE VARINT(X,Y,F,ALPHA, BETA,AREA) 

GENERATING VARIATIONAL FINITE ELEMENT AREA INTEGRATIONS OF THE 
LINEAR BASIS FUNCTION LAGRANGIAN FOR THE HELMHOLTZ EQUATION. 
THESE ARE RETURNED IN F(3,3). X(3) AND Y(3) ARE THE WAVENUMBER 


NORMALIZED CARTESIAN COORDINATES OF THE TRIANGLE VERTICES. 
X = Ko*¥x, Y = Ko*y ALPHA AND BETA ARE COMPLEX 


C MATERIAL PARAMETERS WITHIN THE ELEMENT. 
C 
C FOR TM INCIDENCE: ALPHA = 1/ur; BETA = er 
c FOR TE INCIDENCE: ALPHA = 1/er; BETA = ur 
e 
C OUTPUTS : F(3,3) - FINITE ELEMENT AREA INTEGRATION 
C AREA - AREA OF A TRIANGLE 
@ 
COMPLEX ALPHA, BETA, B12, F(3,3) 
REAL X(3), Y(3), T(3,3), AREA, DET 
INTEGER L, K 
e 


DET = ABS(X(2)*¥(3) + X(3)*¥(1) + X(1)*Y(2) - X(3)*Y(2) - 
CX(1)*¥(3) - X(2)*¥(1)) 

AREA = ABS(0. 5*DET) 

B12 = BETA/12. 


TG) SC eee) DER 

TG) = Cul ea) Daw 

TO CiGh ie tC2) Dak 

pCO) CRS) = 38C2)))) /0lant 

W223) SIRO OK S))))) 11 aT 

WO) SION = ED eda 

13 RSC NC2)) Gs C2 
BMC) ONCE IG) poe ena sy) Noles, 
TS 3 a 2) 2) ee 
DO 10, K = 


80) il} hy = al, 
F(K,L) = ALPHA*(T(1,K)*T(1,L) + T(2,K)*T(2,L)) - B12 
IF(K. EQsn) WCK,L) = FCK.E) = Bie 
10 F(K,L) = AREA*F(K,L) 


C 
RETURN 
END 

C 

C 


Ce a eG? GCC? Ca Cl] Co Cen Cl ce cy Co) OC) Cl Cy CO aI Ia aI Ooo 


APPENDIX E. FIELD FEEDBACK PROGRAM 


FIELD FEEDBACK FORMULATION PROGRAM 
WRITTEN BY T.B. WELCH 


w/ PROGRAMMING IDEAS FROM PROF M.A. MORGAN 


BCOND : 
OFFSET z 
ANS = 
DPER - 


PERND = 


it - 


CNVEC = 
MRES = 


BOUNDARY CONDITIONS 
OFFSET IN WAVELENGTHS (PERIMETER TO BOUNDARY) 
CALCULATED PSI VALUES ON PERIMETER 
DESIRED PERCENT ERROR SCALE FACTOR FOR GREEN'S 
FUNCTION INTEGRAL PATCH STEPPING 
BIAS THAT IS ADDED TO THE GFI STEP FOR < 1.0 
Plas (HAT 1s ADEE TOstHE GF STEP FOR > 1.0 
MAXIMUM DISTANCE BEYOND WHICH NO CONTRIBUTION IS MADE 
TOP Tee cra 
WAVENUMBER 
X ORIGIN 
Y ORIGIN 
NUMBER OF EVENLY SPACED POINTS DESIRED FOR THE FAR 
FIELD CALCULATIONS (360/NRES = ANGULAR RESOLUTION) 
MATRIX THAT RELATES EACH BOUNDARY NODE VALUE TO THE 
UNKNOWN PERIMETER NODE VALUE. MULTIPLY U BY A DRIVING 
VECTOR (ON BOUNDARY) TO FIND PERIMETER VALUES. 
NUMBER OF PERIMETER NODES 
GEOMETRY ARRAY CONTAINING: 
POSITION OF PERIMETER NODES 
POSITION OF PERIMETER NODES 
UNIT NORMAL OF PERIMETER NODES 
UNIT NORMAL OF PERIMETER NODES 
POSITION OF GEOMETRIC CONTOUR NODES 
POSITION OF GEOMETRIC CONTOUR NODES 
POSITION OF BOUNDARY NODES 

Sir = POSITION OF BOUNDARY NODES 
MATRIX THAT RELATES PERIMETER VALUES BACK OUT 
TO THE BOUNDARY VIA A GREEN'S FUNCTION INTEGRAL 
VECTOR OF SCATTERED FIELD BACK ONTO THE BOUNDARY 
MESH RESOLUTION 


Hn & W 
t 
KS KS KEG DM 


COMPLEX BCOND( 100) ,ANS( 100) ,U( 100,100) ,T( 100,100) ,;CNVEC( 100) 
REAL OFFSET ,MESH(100,8) ,DPER,K,XORG, YORG, MRES , MAXD 
INTEGER PERND,I,J,NRES,LBIAS,GBIAS 


OPEN (UNIT 
OPEN (UNIT 


40, FILE 
50, FILE 


'C: MSFORT F3. DAT' ,STATUS=' UNKNOWN’ ) 
'C: MSFORT FFPAT. DAT' ,STATUS=' UNKNOWN’ ) 


CALL INPUT(BCOND, ANS ,U,OFFSET, PERND ,MESH,DPER,K,NRES, XORG, YORG, 
CMRES , LBIAS,GBIAS ,MAXD) 

CALL TMAT(U,PERND,MESH,T,OFFSET,DPER,BCOND ,MRES, LBIAS ,GBIAS , MAXD) 
CALL CNSOLV(T,BCOND,CNVEC , PERND) 


Cae) 


C22 2 (a ©) 


CALL FFLD(CNVEC , PERND,MESH,U,OFFSET,K,NRES, XORG, YORG) 


CLOSE(40) 
CLOSE(50) 


STOP 
END 


SUBROUTINE INPUT( BCOND,ANS ,U,OFFSET , PERND ,MESH, DPER,K,NRES, XORG, 


CYORG ,MRES, LBIAS ,GBIAS ,MAXD) 


THIS SUBROUTINE READS THE FINITE ELEMENT PROBLEM DATA FROM 


THE DATA FILE CALLED 


PRS. DAT sae 


SOLVE THE FIELD FEEDBACK FORMULATION, (F3). 


COMPLEX BCOND( 100) ,ANS( 100) ,U(C 100,100) 


REAL OFFSET,MESH(100,8),DPER,K,XORG, YORG ,MRES , MAXD 


INTEGER PERND,I,J,NRES,LBIAS,GBIAS 


WRITE(* ,**) 


READ(40,*) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,**) 
READ(40,*) 


WRITE(6,*) 
WRITE(6,**) 
WRITE(6,**) 
WRITE(6,**) 
WRITE(6,*) 
WRITE(6,**) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,**) 
WRITE(6,*) 


DOT 10 Sie 


DO 10, 
READ(40,*) MESH(J,I) 


CONTINUE 


DO 20, J = 


" READING INPUT DATA ' 


PERND 
OF Psi] 
DPER 
K 
XORG 
YORG 
NRES 
MRES 
LBIAS 
GBIAS 
MAXD 


' NUMBER OF PERIMETER NODES 
"BOUNDARY CONTOUR OFFSET 

' DESIRED GFI SCALE FACTOR 
' GFI STEP BIAS FOR < 1.0 
' GFI STEP BIAS FOR > 1.0 
' MAX DIST > NO CONTRIBUTION TO GFI 
' WAVENUMBER 
' X ORIGIN 

' Y ORIGIN 

' NUMBER OF NODES FOR SIGMA 
' REQUESTED MESH RESOLUTION 


iL oh 
J = 1, PERND 


1 Pee 


READ(40,*) BCOND(J) 


CONTINUE 


DO 30, J = 


1, PERND 
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READ(40,*) ANS(J) 
CONTINUE 


DO 40, I = 1, PERND 
DO 40, J = 1, PERND 
READ(40,*) UCI,J) 
CONTINUE 


DO 50, I = 1, PERND 


MESH(I,5) = MESH(I,1) + MESH(1I,3)*OFFSET/2. 0 
MESH(I,6) = MESH(I,2) + MESH(I,4)*OFFSET/2. 0 
MESH(I,7) = MESH(I,1) + MESH(I,3)*OFFSET 
MESH(I,8) = MESH(I,2) + MESH(1,4)*OFFSET 
CONTINUE 
RETURN 


END 


SURBROUTINE TMATCU;PERND,MESH,T,OFFSET,DPER, BCOND,MRES,LBIAS,GBIAS, 
CMAXD) 


THIS SUBROUTINE CALCULATES THE GREEN'S FUNCTION INTEGRAL 
(FOR A SINGLE BASIS FUNCTION BOUNDARY CONDITION) GEOMETRIC 
PERIMETER WITH RESPECT TO EACH OF THE OFFSET BOUNDARY NODES. 


THE INTEGRATION IS REPEATED UNTIL EACH BOUNDARY CONDITION HAS BEEN 


INDIVIDUALLY APPLIED AND INTEGRATED WITH RESPECT TO EACH OFFSET 
BOUNDARY NODE. THESE VALUES ARE RETURNED IN THE " T " MATRIX 
FOR USE IN THE FIELD FEEDBACK FORMULATION. THE MATRIX IS 
ORGANIZED, T m,n. 


COMPLEX U( 100,100) ,T( 100,100) ,PVEC( 100) , PDVEC( 100) 

COMPLEX J,HORP1,HORP2,H1RP1,SUM,TEMP( 100) , BCOND( 100) 

COMPLEX H1RP2,PSI,PSIRP,INTEGRAL,DPSIC 

REAL RMRP,NORM11,NORM12,DOT,DPER ,DISTM ,MAXD 

REAL OFFSET,MESH(100,8) ,DIST,R,DZ,DL,MRES 

INTEGER 1,M,N,NN,STEP,FN,SN,PERND ,MM,STEPMX ,STEPMN, LBIAS ,GBIAS 
OPEN (UNIT = 2, FILE = 'C: MSFORT TMAT. DAT’ ,STATUS = 'UNKNOWN' ) 


J = (0.0,1.0) 


STEPMX = INT(DPER*(-48. 0*OFFSET + 17.2) + LBIAS) 

STEPMN = INT(DPER + GBIAS) 

WRITE(*,*) ' LOADING T MATRIX ' 

WRITE(*,*) ' MAXIMUM STEP = ',STEPMX,', MINIMUM STEP = ' ,STEPMN 
WRITE(*,*) ' MAXIMUM DISTANCE FOR ANY CONTRIBUTION = ‘ ,MAXD 

DO 40, M = 1, PERND 


DO 5, I = 1, PERND 
IF(M.EQ. 1) THEN 
PVEC(I) = (1.0 + UCI,M))/2.0 
PDVEC(I) = (1.0 - U(I,M))/OFFSET 
ELSE 
PVEC(I) = U(I,M)/2.0 
PDVEC(I) = -U(I,M)/OFFSET 


ENDIF 
CONTINUE 


DO 30, N = 1, PERND 
WRITE(*,1000) M, N, PERND 
SUM = (0.0,0.0) 

DO 20, NN = 1, PERND 
FN = NN 
IF(NN. EQ. PERND) THEN 
SN = 1 
ELSE 
SN = NN+1 
ENDIF 
DIST = SQRT((MESH(FN,5)-MESH(SN,5) )**2+(MESH(FN, 6) 
-MESH(SN, 6) )**2) 
INTEGRAL = (0. 0,0.0) 


DISTM = SQRT( CMESH(N, 7) -MESHCFN,5) )**2+( MESH(N, 8) =MESH( EN, Gi e2s 


R = MESH(SN,5) - MESH(FN,5) 
DZ = MESH(SN,6) - MESH(FN,6) 
DL = SQRT(R**2 + DZ**2) 
NORM11 = -DZ/DL 

NORM12 = R/DL 


IF(CDISTM. GT. MAXD) THEN 
GOT y 2G 

ELSE IE CDIS PS EE lope CHEN 
STEP = STEPMA 

ELSE 
STEP = STEPMN 

ENDIF 

IFCSTEP. LI. 1) estee sed 


DO 10, I = 1, STEP+1 

IF(I.EQ.1) THEN 
RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+0. 25**(MESH(SN,5)- 
MESH(FN,5))/FLOAT( STEP) ) )***2+(MESH(N,8)-(MESH(FN,6)+ 
0. 25°*(MESH(SN,6)-MESH(EN, 6) ) /FLOAT( STEP) ) )***2) 
DPSIC = PDVEC(FN) + 0. 25*( PDVEC(SN) -PDVEC(FN) )/STEP 
PSIRP = PVEC(FN) + 0. 25%(PVEC(SN)-PVEC(FN) )/STEP 
DOT = (NORM11**(MESH(N,7) - (MESH(FN,5)+0. 25**(MESH(SN,5)- 
MESH(FN,5))/STEP) )+(NORM12**( MESH(N, 8) -(MESH(FN,6)+0. 25% 
(MESH(SN,6) -MESH(FN,6))/STEP) )))/RMRP 

ELSEIF(I. EQ. STEP+1) THEN 
RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+( FLOAT( STEP) -0. 25)*( 


MESH(SN,5)-MESH(FN,5) )/FLOAT( STEP) ) )***2+( MESH(N, 8) -(MESH 
(FN ,6)+(FLOAT( STEP) -0. 25)*(MESH(SN,6)-MESH(FN, 6) )/FLOAT 
(STEP) ) )7*2) 

DPSIC=PDVEC( FN)+( FLOAT( STEP) -0. 25)*( PDVEC( SN) -PDVEC( FN) ) 
/(STEP) 

PSIRP = PVEC(FN)+( FLOAT( STEP) -0. 25)*( PVEC(SN) -PVEC(FN) ) / 
(STEP) 


DOT = (NORM11*( MESH(N, 7) -(MESH(C EN, 5 )+(@EOAT(STEP)—OGn2Zay5. 
(MESHCSN,5) -MESH( EN, 5)')/STEP) )-( NORMAZ* ( MESH(N 76> @iiem 
(FN,6)+( FLOAT(STEP)-0. 25)*"(MESH(SN,6)-MESH( FN, 6)))/> Ueno 
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))) /RMRP 

ELSE 
RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+FLOAT( I-1)*(MESH(SN, 
5) -MESH(FN,5))/FLOAT( STEP) ) )***2+(MESH(N,8)-(MESH(FN,6)+ 
FLOAT( I-1)**(MESH(SN,6) -MESH(FN,6))/FLOAT( STEP) ) )****2) 
DPS IC=PDVEC( FN)+( I-1)**( PDVEC( SN) -PDVEC( FN) )/( STEP) 
PSIRP = PVEC(FN)+(I-1)*( PVEC(SN) -PVEC(FN) ) /( STEP) 
DOT = (NORM11**(MESH(N,7)-(MESH(FN,5)+(I-1)*(MESH(SN,5)- 
MESH( FN,5))/STEP) )+(NORM12**( MESH(N,8)-(MESH(FN,6)+(I-1)* 
(MESH(SN,6)-MESH(FN,6))/STEP) )))/RMRP 

ENDIF 

CALL HAN1(CMPLX(RMRP,0.0),HORP1,H1RP1) 

IF((I. EQ. 1). OR. (I. EQ. STEP+1)) THEN 
PSI = (J/4.0)*(HORP1*DPSIC - PSIRP*DOT*H1IRP1)/2 


ELSE 
PSI = (J/4.0)*(HORPI*=DPSIC = PSIRP*DOT*HIRP1) 
ENDIF 
INTEGRAL = INTEGRAL + PSI*DIST/STEP 
CONTINUE 
SUM = SUM + INTEGRAL 
CONTINUE 
TCN,M) = SUM 
CONTINUE 
CONTINUE 


DO 55, I = 1, PERND 
TEMP(I) = (0.0,0.0) 
DO 50, M = 1, PERND 
TEMP(I) = TEMP(I) + T(1I,M)**BCOND(M) 
CONTINUE 
WRITE(2,%) (TC(I,MM), MM = 1, PERND) 
CONTINUE 


DO 60, I = 1, PERND 
WRITE(2,*) I, TEMP(I) 
CONTINUE 
AOrhMolGix me COLUMN) '.13, , ROW ',13,;' OUT OF ‘,I3) 
CLOSE(2) 
RETURN 
END 
SUBROUTINE CNSOLV(T,BCOND , CNVEC, PERND) 
THIS SUBROUTINE CALCULATES THE C VECTOR BY SOLVING: 


= 
Cn = [I - T] * BOUNDARY CONDITIONS (INCIDENT FIELDS) 


COMPLEX BCOND(100),T(100,100),TEMP( 100) ,CNVEC( 100), DETERM 
REAL COND, DMAG 


ea 


INTEGER PERND {1 {JK ib, 1NCKiAaNiEe. 


INORM = 0 
NMAX = 100 


DO 10, I = 1, PERND 
DO 10, J = 1, PERND 
IF(I. EQ. J) THEN 
T(1, Jo = 1 =e) 
ELSE 
TCIgd) = -TGie) 
ENDIF 
10 CONTINUE 
WRITE(*,*) "INVERTING THE [I - T] MATRIX 


@ 
CALL CSMINV(T,NMAX, PERND , DETERM ,COND , INORM) 
C 
DMAG = CABS(DETERM) 
WRITE(*,*) ' DETERMINANT = ', DMAG 
WRITE(*,**) ' CONDITION NUMBER = ', COND 
WRITE(*,*) ' MULTIPLING MATRICES TO FORM THE SCATTERED FIELDS ' 
C 


DO=30 7) 1 =v EeRNG 
CNVECCT )" =" ( 020.020) 
DON30, J = 12S eERND 
CNVEC(I) = CNVEC(I) + TC I,J)*BCOND( J) 
30 CONTINUE 


C 
RETURN 
END 
C 
G 
SUBROUTINE FFLD(CNVEC, PERND ,MESH,U,OFFSET,K,NRES , XORG, YORG) 
: \ 
CG THIS SUBROUTINE CALCULATES THE FAR FIELDS DUE TO THE OFFSET 
C BOUNDARY SCATTERED FIELDS AND THE PERIMETER SCATTERED FIELDS. 
G ADDITIONAL GREEN'S FUNCTION INTEGRALS ARE ACCOMPLISHED. 
@ 
COMPLEX CNVEC( 100) ,U(100,100),J,PSISP( 100) ,TEMP,PSI,DPSI 
COMPLEX DPSISP( 100) , INTEGRAL,DPSIM( 100,100) ,PSIM( 100,100) 
REAL MESH(100,8),OFFSET,K,DOT,DOT1,PI,ARES,XORG, YORG,DIST,SIGMA 
REAL R,DZ,DL,NORM11,NORM12 
INTEGER PERND,I,M,N,L,NRES,FN,SN,STEP,II 
C 
WRITE(*,*) ' CALCULATING THE SCATTERED FIELDS’ 
J =5(0705 be) 
PI = 4. O*ATAN(1. 0) 
ARES = 2. 0*PI/FLOAT(NRES) 
STEP = 5 
‘@ 
DO 5, M = 1, PERND 
DO 1, PERND 


PSIM(I,M) = (1.0 + UCI,M))/2.0 
DPSIM(I,M) = (1.0 - U(1I,M))/OFFSET 
ELSE 
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PSIM(I,M) = U(I,M)/2.0 
DPSIM(I,M) = -U(1I,M)/OFFSET 
ENDIF 
CONTINUE 


DO 10, I = 1, PERND 
PSISP(I) = (0.0,0.0) 
DPSISP(1) = (0.0,0.0) 
DO 10, L = 1, PERND 
PSISP(I) = PSISP(I) + PSIM(I,L)*CNVEC(L) 
DPSISP(I) = DPSISP(I) + DPSIM(1,L)*CNVEC(L) 
CONTINUE 


DO 30, I = 0, NRES-1 
WRITE(6,1000) I+1, NRES 
INTEGRAL = (0.0,0.0) 

DO 20, N= 1, PERND 


FN = N 
IF(N. EQ. PERND) THEN 
SN = 1 
ELSE 
SN=N+1 
ENDIF 


R = MESH(SN,5) - MESH(FN,5) 
DZ = MESH(SN,6) - MESH(FN,6) 
DL = SQRT(R**2 + DZ*2) 
NORM11 = -DZ/DL 
NORM12 = R/DL 
DO 15, II = 1, STEP+1 
DIST = SQRT((MESH(FN,5)-MESH(SN,5) )**2+( MESH(FN,6)- 
MESH(SN, 6) )****2) 
IF(II.EQ.1) THEN 
PSI = PSISP(FN)+0. 25*( PSISP(SN) -PSISP( FN) )/ 
FLOAT( STEP) 
DPSI = DPSISP(FN)+0. 25*(DPSISP( SN) -DPSISP 
(FN) )/FLOAT( STEP) 
DOT = NORM11*SIN(I*ARES) + NORM12*COS(I**ARES) 
DOT1 = (MESH(FN,5)+0. 25**(MESH(SN,5)-MESH(FN, 
5))/FLOAT( STEP) -XORG)*SIN(I*ARES) + (MESH(FN,6)+ 
0. 25**(MESH(SN,6)-MESH(FN,6) ) /FLOAT( STEP) - 
YORG)**COS( I**ARES) 
ELSEIF(II. EQ. STEP+1) THEN 
PSI = PSISP(FN)+(FLOAT( STEP) -0. 25)*( PSISP(SN) - 
PSISP(FN) ) /FLOAT( STEP) 
DPSI = DPSISP(FN)+( FLOAT( STEP) -0. 25)*(DPSISP(SN)- 
DPSISP(FN) ) /FLOAT( STEP) 
DOT = NORM11*SIN(I*ARES) + NORM12*COS( I*ARES) 
DOT1 = (MESH(FN,5)+(FLOAT( STEP) -0. 25)*(MESH(SN,5)- 


MESH(FN,5))/FLOAT( STEP) -XORG)*SIN(I**ARES) + (MESH( 
FN, 6)+(FLOAT( STEP) -0. 25)**(MESH( SN, 6) -MESH(FN,6))/ 
FLOAT( STEP) -YORG)*COS(I**ARES) 

ELSE 
PSI = PSISP(FN)+FLOAT( II-1)*( PSISP(SN)-PSISP(FN) ) / 
FLOAT( STEP) 
DPSI = DPSISP(FN)+FLOAT( II-1)*(DPSISP( SN) -DPSISP 
(FN) ) /FLOAT( STEP) 


DOT = NORM11*SINCI*ARES) + NORM12*COS( I*ARES) 

DOT1 = (MESH(FN ,5)+FLOAT( II-1)*(MESH(SN,5)-MESH(FN, 
5))/FLOATC STEP) -XORG)*SIN(1* ARES Je (MESH ENecre. 
FLOAT( II-1)*(MESH(SN, 6) -MESHCFN,6))/FLOAT( STEP) - 


13) 
20 
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maa 
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YORG)**COS( I*ARES) 
ENDIF 


IFC( 11. EQ. 1). OR. C11. EQsS1i Pe) meet 


TEMP=( J*DPSI+DOT*PSI )*(EXP( J*DOT1) )*DIST/(2. 0* 


C STEP) 
ELSE 


TEMP=(J*DPSI+DOT*PS? )*(EXP(J*D0n1) ) bisa, 


C ColBE 
ENDIF 
INTEGRAL = INTEGRAL + TEMP 
CONTINUE 
CONTINUE 
SIGMA = ((CABS( INTEGRAL) )**2.0)/(4. O*K) 


WRITE(50,1010) I+1, (I *ARES*180. 0/7P1), SIGMA 


CONTINUE 


WRITE(6,*) ' 


WRITE(6,*) ' <<< FIELD PATTERN STORED IN FFPAT. DAT >>> ' 


WRITE(6,*) ' ' 


FORMAT( 1X, ‘INTEGRAL "13, , OUD OF 13) seGtPhETED.) 


FORMAT(1X%,13,2x%,F6. 2 22 aes) 


RETURN 
END 


SUBROUTINE HAN1(Z,HO,H1) 


Computing Hankel Functions for n=0,1 with 


Complex Argument, Z. Direct Power Series Method for 
CABS(Z) .LE. 5 and Hankel's Asymptotic Formula for 
CABS(CZ)n.. Gil on Written 11/6/87 by M.A. Morgan 


INTEGER M,M2 
REAL C(34) ,DM,F( 34) ,GO,P(34) ,Pi, P2 


COMPLEX Z,Z2,Z3,24,J0,J1,Y0,Y1,AM,CL, PO,P1,00,Q1 


COMPLEX EO,E1,X0,X1,HO,H1, j 
PI=3. 1415927 

P2=2.0/PI 

7200.1) 
IF(CABS(Z). LE. 5.0) THEN 


Direct Power Series Method 
GO= 1.78072 
Z2=0. 5*Z 
CL=CLOG( G0*Z2 ) 
Computing F(m) = m! and Pin) = 1 tee ea 


F(1l)=i96 


130 


ieee te mM 


LSD 1 | 


P(1)=1.0 

DO 11 M=2,34 
F(M)=M*F(M-1) 
P(M)=P(M-1)+1. 0/M 

CONTINUE 


Computing Power Series Coefficients 


DM=-1. 0 

DO 22 M=1,34 
C(M)=DM/(FCM)*F(M) ) 
DM=-DM 

CONTINUE 


Computing JO and Jl 


JO=(1. ,0. ) 

J1=(0. ,0. ) 

M=0 

M=M+1 

M2=27M 

AM=C(M)*( Z27"M2 ) 

JO=JO+AM 

J1=J1-M*AM 

IF((CABS( AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 33 
n= oy 2 


Computing YO and Yl 


M=0 
YO=CL**JO 
Y1=Z2*CL*J1-0. 5*JO 
M=M+1 
M2=2"M 
AM=C(M)**P(M)**( Z2%**M2) 
YO=YO-AM 
Y1=Y1+M"'AM 
IF((CABS(AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 44 
YO=P2**Y0 
Y1=P2"Y1/Z2 
HO=JO- j**YO 
H1=J1-4*Y1 
RETURN 

ELSE 


Hankel’ Asymptotic Formula (Abram. & Stegun p. 364 


7 2= 757 

Z3=Z*Z2 

Z4=Z*Z3 

PO=1. 0-. 0703125/2Z2+. 1121521/24 
QO=-. 125/Z+. 0732422/Z3 

P1=1. 0+. 1171875/Z2-. 1441956/24 
Q1=. 375/Z-. 10253906/Z3 

XO=(Z-. 25*PI) 

xl=(2- 7 oP) 
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EO=CEXP( - j**X0) 

E1=CEXP( - j**X1) 

AM=CSQRT(P2/Z) 

HO=AM**( PO-4**Q0)**E0 

H1=AM*(P1-9*Q1)*E1 
ENDIF 


RETURN 

END 

SUBROUTINE CSMINV(A,NDIM,N,DETERM, COND, INORM) 

INORM - FLAG TO NORMALIZE COLUMNS AND ROWS OF MATRIX A 


MATRIX NORMALIZATION BY M.A. MORGAN 
APRIL 24,1978 


A - MATRIX TO INVERT = INPUT /OUTE Un 
NDIM = UN POT 

N - = NP Or 

DEEKM = DE TERMINATE OF TA =" OUDEUR 

COND - CONDITION NUMBER OF A = OUTrer 


INORM - INTEGER NORMALIZATION FLAG - INPUT 


COMPLEX A( 100,100) ,PIVOT( 100) ,AMAX,T,SWAP , DETERM, U 

INTEGER 1,J,K,L,1PIVOT( 100) , INDEX( 100, 2), [ROW 1COGLUM ie ay 
INTEGER JCOLUM ,N, INORM 

REAL TEMP ,ALPHA( 100) ,COL( 100) ,ROW(C 100) ,AJK, SUMAXA, SUMROW , SUMAXI 


IF(NDIM. GT. 100) THEN 
WRITE(*,*) ' ERROR IN INVERTION CALL... DIMENSION > 100 ' 
STOP 

ENDIF 

IF(N. GT. NDIM) THEN 
WRITE(*,*) ' ERROR IN INVERTION CALL... N > MAX DIM. ' 
STOP 

ENDIF 

IF( INORM. Nie) GO sone 


0.0 
1,N 
AJK = CABS(A(J,K)) 
IF( AJK. GT. COL(K)) COL(K) = AJK 
CONTINUE 
Der 2 — see 
A(J,K) = ACJ,K)/COL(K) 
CONTINUE 
ROW NORMALIZING 
pO6J=1,N 
ROW(J) = 0.0 
DO 4K = 1, N 
AJK = CABS(A(J,K)) 
IF(AJK. GT. ROW(J)) ROW(J) = AJK 
CONTINUE 


=, 
O 
Oo 
A 
TL Zt 
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200 


260 


DO5 K=1, 
ACN) 


AC 

CONTINUE 

CONTINUE 

DEGERM = CMPLA(C1.0,0. 

SUMAXA = 0.0 

BO 20 J =1, N 
ALPHA(J) = 0.0 
SUMROW = 0. 
DO 10 I = _ N 

ALPHA( J) = 

SUMROW = SUMROW 
ALPHA( J) = SQRT( 


J,K)/ROWC J) 


0) 


ALPHA(J) + ACJ,1)* CONJGCACJ,I)) 
+ CABS(CACJ,I1)) 
ALPHA(J)) 


IF(SUMROW. GT. SUMAXA) SUMAXA = SUMROW 


IPIVOT(J) = 0 
DO 600 I= 1, N 


AMAX = CMPLX(0. 0 
be 105 J ="1, N 
IF (IPIVOT( 

DO 100 K = 

ie (Give 


BO 0) 


JJ OO. 10560 
eal 
TVOT(R-1) 80; 100, 740 


TEMP = AMAX*CONJG(C AMAX) - AC J,K)*CONJGCACJ,K)) 


IF (TEM 


CONTINUE 
CONTINUE 


eee 385,100 


ier OUCiTCOLUM). = IPIVOTCICOLUM) + 1 


IF (IROW-ICOLUM) 
DETERM = -D 
NOe200 Le = 

SWAP = 


ACIROW, 


AC ICOLUM,1) 
SWAP = ALPH 
ALPHA( IROW) 
ALPHA( ICOLU 
INDEX( 1,1) 


INDiExXC I, 2) >= 


PIVOT(I) = 
U = PIVOT(I 
DETERM = DE 
DETERM = DE 


140, 260, 140 
ETERM 

1, N 
AC IROW,L) 
L) = ACICOLUM,L) 
= SWAP 
A( IROW) 

= ALPHA( ICOLUM) 
M) = SWAP 
= IROW 
ICOLUM 
A( ICOLUM, ICOLUM) 
) 
TERM*U 
TERM/ALPHA( ICOLUM) 


TEMP = PIVOT(1)*CONJG( PIVOT(I1)) 


IF( TEMP) 330, 


A( ICOLUM, IC 
10) S30 1h = 

= Pi 
AC ICOLUM, L) 


720750950 


OLUM) = CMPLX(1. 0,0. 0) 
1, N 

VOT(1) 

= ACICOLUM,L)/U 


los 


OS 
710 


900 


210 


DO 550 Li = 1, N 
IF(L1-ICOLUM) 400, 550, 400 
T = A(L1,ICOLUM) 
A(L1,ICOLUM) = CMPLX(0. 0,0. 0) 
DO 450 L = 1, N 
U = AC ICOLUM,L) 
ACA) = AC en 


CONTINUE 
CONTINUE 
DO 710 T= 1, N 
L = Nie ee 
TF CINDEXCES 1!) —INDEXGI 2G 20 7 tOnmcs! 


JROW = INDEX(L,1) 
JCOLUM = INDEX(L, 2) 
DO? 05 eke — eee 
SWAP = A(K, JROW) 
ACK, JROW)s= ACK ICOLU 
AC K,JCOLUM) = SWAP 
CONTINUE 
CONTINUE 
SUMAXI = 0.0 
DO 910 I =1, N 
SUMROW = 0.0 
DO 900s) —— ian 
SUMROW = SUMROW + CABS(A(I,J)) 
IF(SUMROW. GT. SUMAXI) SUMAXI = SUMROW 
CONTINUE 
COND = SUMAXA*SUMAXI 
IFC INORM. NE.1) GO TO 955 
DO 950 K = 1, N 
HO 9500 = ea 
ACJ,K) = ACJ,K)/CROW(K)*COL(J)) 
CONTINUE 
RETURN 
WRITE(*, 730) 
FORMAT( ' MATRIX IS SINGULAR’ ) 
RETURN 
END 
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APPENDIX F. DIELECTRIC CYLINDER SCATTERING PROGRAM 
ccccccccccccccecccccccccccccccecccecccccccccccccccccccccccccccecccccccce 


C 
PLANEWAVE SCATTERING BY A DIELECTRIC CYLINDER C 
E - WAVE (TM CASE) C 
H - WAVE (TE CASE) C 
C 
C 


Circ? Cre 


cCcccccccccccccccccccCcccccccccccccccccccccccccccccccccccccccccccccccccc 


© 


COMPLEX GAMMA(0: 200,2),SIGMA(2),ER,MU,KR,YR,ZR 
COMPLEX*16 JA(0: 200) ,DJA(0: 200) ,KRRA 

COMPLEX*16 J(0: 200) ,DJ(0: 200) 

REAL**8 Y(0: 200) ,DY(0: 200) , YA(0: 200) ,DYA(0: 200) , JB(0: 200) 
REAL**8 DJB(0: 200) ,RA,KORA 

REAL KO,PI,SIGMAN(2),A,B 

INTEGER I,1I,MODE,ARES,N 


OPEN(3,FILE='C: MSFORT DECTEM. DAT' ) 
Ee— 3, 1415927 


WRITE(*,*) ‘PERMITTIVITY FORMAT IS "“ a+ jb, " ' 
WRITE(*,**) ‘Enter Dielectric Constant (REAL PART, a) ' 
READ(*,**) A 

WRITE(*,**) ‘Enter Dielectric Constant (IMAGINARY PART, b) ' 
READ(*,**) B 

ER = CMPLX(A,B) 


WRITE(*,*) "PERMEABILITY FORMAT IS " a+ jb, " ' 

WRITE(*,*) ‘Enter Permeability Constant (REAL PART, a) ' 
READ(*,**) A 

WRITE(*,*) ‘Enter Permeability Constant (IMAGINARY PART, b) ' 
READ(*,*) B 

MU = CMPLX(A,B) 


WRITE(*,*) ‘Enter the wave number (ko) ' 
READ(* ,**) KO 


WRITE(*,*) ‘Enter Cylinder Radius (IN WAVENUMBER UNITS)  ' 
WRITE(*,*) ‘WARNING: Do Not Enter Zero ! 
READ(**,**) KORA 


CSQRT(MU*ER) 
CSQRT(ER/MU) 
CSQRT(MU/ER) 
KORA/KO 
KRRA = KR*RA 


K< 
wy 
hou ue 


WRITE(*,*) 'Enter No. of Modes: 
READ(**,**) MODE 

WRITE(*,**) ‘Enter the angular resolution 
READ(*,*) ARES 


es 


WRITE(3,*) ‘Cylinder Scattering vs. angle' 
WRITE(3,110) ER,MU,RA,MODE,KO,KRRA 
WRITE(3,*) 


CALL BES(MODE,KORA,JB,Y,DJB,DY) 
CALL DCBJNS (DCMPLX(KORA,0. OD+00) ,MODE,J,DJ) 
CALL DCBJNS ((KRRA*KO),MODE,JA,DJA) 
WRITE(**,**) ' RETURNED FROM FINAL BESSEL CALL’ 


CALCULATING GAMMAs 


Ca-Ga.C9 


DO 10, N= 0, MODE 
GAMMA(N,1) = (JAC(N)*DJ(N)-YR*DJACN)*J(N) )/( JACN)*CMPLX(DJ(N) 
C ,-DY(N)) - YR*DJA(N)*CMPLX(J(N) ,-Y(N))) 
GAMMA(N,2) = (JA(N)*DJ(N)-ZR*DJA(N)*J(N) )/( JACN)*CMPLX(DJ(N) 
6 ,-DY(N)) - ZR*DJA(N)*CMPLX(J(N),-Y(N))) 
WRITE(**, 1000) N,GAMMA(N,1),GAMMA(N, 2) 
0 CONTINUE 


CALCULATING SIGMAs 


22 Oe 


DO 30, I = 180, -180, -ARES 
DO 40) I) = ieee 
SIGMA(II) = (0.0,0.0) 
DO 20, N = 1, MODE 
SIGMA(II) = SIGMA(II)+2. O*GAMMA(N, II)*COS(N*PI*I/180. 0) 
20 CONTINUE 
SIGMA(II) = SIGMA(II) + GAMMA(0O,IT) 
SIGMAN(II) = ((4.0/KO)**(CABS( SIGMA(II)))**2) 
40 CONTINUE 
WRITE(3,120) I, SIGMAN(1), SIGMAN(2) 
30 CONTINUE 


6 
6 
110 FORMAT(1X,'Er = ih 6s Gn be Romane 
Co Ma = | FO. Geno. Gy 5 
C' RADIUS (METERs) = fF Bass 
C' MAX MODE = gelcreie. 
G7KO = USES See 
Ciera = ' |F8.5,1X,F8. 5) 


120 FORMAT(1X,14,2(3X,E14.4)) 
1000 FORMAT(1X,13,1X,2(E12.4,1X,E12. 4,4X)) 
e 
C 
STOP 
END 
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Integrated Software Systems Corporation 
10505 Sorrento Valley Road 
San Diego, CA 92121 


meeKVE-DIGITIZER 


West Coast Consultants 
4202 Genesee Avenue, Suite 309 
san Diego, CA 92117 


merviicrosoit FORTRAN 


16011 NE 36th Wavy 
BOX 97017 
Redmond, WA 98073 


. Microwav NDP FORTRAN 
POB 79 
Kingston, MA 023064 


maro!. vl.A. Viorgan, Code 62Mw 
Naval Postgraduate School 
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